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Foreword 


The  relationship  between  sun  glitter  in. the  sea  surface  and  the  statistical 
distribution  of  sea  surface  slopes  allows  oceanographers  to  gather  informa¬ 
tion  on  the  physical  conditions  of  the  sea  surface  by  analyzing  the  glitter. 
Aerial  photographs  have  been  a  convenient  means  for  recording  the  glitter 
pattern,  but  analysis  of  the  photographs  is  often  labor-intensive. 

Current  digital  image  processing  systems  and  techniques  have  reduced  the 
labor  factor  considerably  and  have  increased  practicality  of  sun  glitter  analysis. 
This  report  describes  the  approach  that  was  used  by  the  Naval  Ocean  Research 
and  Development  Activity,  Remote  Sensing  Branch,  to  perform  sun  glitter 
analysis. 


R.  P.  Onorati,  Captain,  USN 
Commanding  Officer,  NORDA 


Executive  summary 


The  image  of  the  sun  reflected  from  the  rough  surface  of  the  sea  forms 
a  diffuse  pattern,  whose  details  depend  on  the  nature  of  the  surface  waves 
and  swell.  Consequently,  statistics  of  the  slope  distribution  of  the  sea  surface 
are  related  to  the  statistics  of  the  sun’s  glitter  on  the  sea  surface.  Aerial 
photographs  are  a  convenient  medium  for  recording  the  glitter  pattern.  The 
mathematical  relationship  to  derive  the  sea  surface  slope  statistics  can  be  deter¬ 
mined  from  an  analysis  of  the  imaging  geometry.  However,  analysis  of  the 
photographs  can  be  a  labor-intensive  procedure. 

The  problem  was  first  studied  over  thirty  years  ago.  Now,  there  are  modem 
digital  image  processing  systems  and  techniques  that  greatly  increase  the  prac¬ 
ticality  of  the  analysis.  This  report  derives  the  relevant  equations  and  describes 
an  implementation  on  the  Interactive  Digital  Satellite  Image  Processing  System 
(IDSIPS).  The  IDSIPS  system  is  operated  by  the  Remote  Sensing  Branch  of 
the  Naval  Ocean  Research  and  Development  Activity  (NORDA).  The  report 
includes  full  formal  documentation  of  the  computer  software.  -*  - 

NORDA ’s  Remote  Sensing  Branch  uses  image  processing  methods  to  ob¬ 
tain  quantitative  information  on  oceanographic  parameters  from  analyses  of 
remotely  sensed  data.  This  work  is  directed  toward  satisfaction  of  U.S.  Navy 
requirements  for  accurate  oceanographic  information. 
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Ocean  wave  slope  statistics  from  automated 
analysis  of  sun  glitter  photographs 


1.  Introduction 

The  statistical  distribution  of  sea-surface  slopes  is  related 
to  the  statistics  of  the  sun’s  glitter  on  the  sea  surface, 
so  information  on  the  physical  conditions  at  the  sea  sur¬ 
face  can  be  obtained  by  analyzing  the  glitter.  Aerial  photo¬ 
graphs  are  a  convenient  medium  for  recording  the  glitter 
pattern.  The  problem  was  studied  several  years  ago  by 
Cox  and  Munk  [1,  2];  their  analysis  is  included  in 
Kinsman’s  book  [3]. 

Modern  digital  image  processing  equipment  and  tech¬ 
niques  afford  a  method  of  performing  the  analysis  in  a 
way  that  is  much  less  labor-intensive.  Another  applica¬ 
tion  of  digital  image  processing  to  photographic  photome¬ 
try  is  reported  in  Reference  4.  This  report  describes  the 
approach  that  was  used  by  the  Naval  Ocean  Research  and 
Development  Activity  (NORDA)  Remote  Sensing  Branch 
to  perform  sun  glitter  analysis.  A  detailed  development 
of  the  formulation  is  followed  by  formal  documentation 
of  the  computer  program. 

2.  Derivations  of  equations 

Geometric  considerations  relate  each  point  in  the  glit¬ 
ter  image  to  the  particular  sea-surface  dlt  and  orientation 
that  are  required  to  reflect  the  sun’s  rays  toward  that  im¬ 
age  point.  The  brightness  of  a  picture  element  (pixel)  is 
related  to  the  fraction  of  the  sea-surface  area  in  the  held 
of  view  that  satisfies  the  geometric  conditions  for  imag¬ 
ing  to  that  pixel.  For  a  given  position  of  the  sun,  the  pix¬ 
el  value  (brightness)  distribution  is  determined  by  the  sea 
conditions.  The  Cox  and  Munk  analysis  infers  informa¬ 
tion  about  sea  conditions  from  analysis  of  the  pixel  value 
distribution  [1-3]. 

2.1  Geometry 

Consider  a  right-hand  Cartesian  coordinate  system 
X1-X2-X3  with  origin  at  the  sea  surface,  X3  vertically  up¬ 
ward,  and  X2  horizontal  and  toward  the  sun  (Refs.  2  and 
3  and  Fig.  1  illustrate  this  geometry).  An  incident  ray 
from  the  sun  is  reflected  from  the  origin  to  a  point  on 
a  horizontal  photographic  plate  above  the  sea  surface.  (The 
case  of  a  photographic  plate  that  is  not  horizontal  is  con¬ 
sidered  in  Section  2.9.)  An  “image”  coordinate  system 


x  y  in  the  plane  of  the  photographic  plate  has  the  >-axis 
pointing  away  from  the  sun.  The  origin  of  the  image 
system  is  the  intersection  of  a  vertical  through  the  center 
of  the  camera’s  aperture  with  the  photographic  plate. 

Define  the  following  angles: 

a  =  direction  of  steepest  ascent  of  the  (element  of)  sea 
surface,  measured  clockwise  from  the  sun  direction  X2. 

(3  =  tilt  of  sea  surface. 

<t>  =  sun’s  elevation  (or  altitude)  above  the  horizontal. 

fi  =  angle  between  vertical  and  vector  from  center  of 
camera’s  aperture  to  image  point. 

v  =  angle  in  image  plane  between  y-axis  and  vector 
from  x-y  origin  to  image  point,  positive  clockwise  from 
y-axis. 

u  =  angle  of  incidence  or  reflection. 

By  applying  the  geometrical  optics  law  of  reflection  [5] 
to  the  reflection  by  the  sea  surface  of  a  ray  from  the  sun 
to  the  photographic  plate,  it  is  found  that  [2,  3] 


cosu  =  cosff  sintp  -  cosa  sin  (5  cos<p,  (1) 

cosfi  =  2cos(3  cosu  -  simp,  (2) 

cot  v  =  cot  a  +  'A  esc  a  cscfi  seco)  cos<p.  (3) 

The  angles  fi  and  v  can  easily  be  calculated  for  any 
image  point.  The  solar  elevation  <p  can  be  calculated  from 
external  conditions  by  standard  techniques.  It  will  be  shown 
in  Section  2.4  that  the  angle  of  incidence  u  can  be  deter¬ 
mined  without  reference  to  a  and  jS.  So  the  mathematical 
problem  consists  of  finding  a  and  |8  for  each  image  point 
from  the  other  conditions,  and  constructing  the  histograms 
that  approximate  the  desired  probability  distributions.  It 
will  be  shown  in  Section  2.7  and  Appendix  A  that 
References  2  and  3  have  a  sign  error  in  the  equation  cor¬ 
responding  to  Equation  (3). 

It  is  convenient  to  define  another  coordinate  system  in 
the  image  plane.  The  system  •X/  j'/  has  the  same  origin  as 
xy,  but  the  Jt  ,-axis  points  toward  the  starboard  wing  of 
the  aircraft  and  the  y  ,-axis  points  toward  the  nose.  Con¬ 
sequently,  the  y,  -axis  points  toward  the  top  of  the  im¬ 
age,  with  the  camera  mounting  geometry  that  was  used. 
Figure  2  illustrates  the  two  coordinate  systems.  It  is  seen 
that  the  x  y  system  is  oriented  at  an  angle  (Az  +  Hdg) 
counterclockwise  with  respect  to  the  x,-y,  system,  where 
Az  =  solar  azimuth,  measured  counterclockwise  from  S, 
Hdg  =  aircraft  heading,  measured  clockwise  from  N. 

Figure  2  also  illustrates  an  image  point  P  and  its  associated 
angle  v. 


2.2  Solar  elevation  and  azimuth 

The  first  step  in  the  analysis  is  to  locate  the  sun  in  the 
imaging  geometry.  Let 

t  =  the  time  of  day  (GMT)  past  midnight  in  minutes, 
X  =  the  longitude  west  of  Greenwich  in  degrees. 
The  local  hour  angle  h  of  the  sun  (the  angular  distance 
of  that  body  west  of  the  local  celestial  meridian)  [6]  is, 
in  degrees, 

h  =  [(t  -  7 20)/ 4]  -  X.  (4) 

The  sun’s  elevation  and  azimuth  can  be  found  in  terms 
of  the  hour  angle  by  solution  of  the  navigational  triangle 
(see  Fig.  3).  The  navigational  triangle  is  a  spherical  triangle 
defined  by  three  points  on  the  surface  of  the  earth  and 
formed  by  the  arcs  of  the  great  circles  that  connect  these 
points  [6].  The  three  points  are  the  position  of  the  observer, 
M,  the  geographical  position  of  the  celestial  body  (the  sun 
in  this  case)  being  observed,  GP,  and  the  earth’s  pole 
nearer  the  observer,  P.  Figure  3  shows  the  navigational 
triangle  with  the  sides  and  one  of  the  vertex  angles  labeled. 
In  Figure  3 

6  =  declination  angle  of  the  sun, 

L  =  latitude  of  the  observer. 

Then  the  result  of  applying  one  of  the  standard  formulas 
for  the  solution  of  spherical  triangles  gives  [7] 

sin<P  =  sinb  sin  L  +  cosh  cos  L  cos  h.  (5) 

The  angle  at  M  is  Az  -  180°,  as  azimuth  was  defined 
in  Section  2.1.  So  the  application  of  another  spherical 
trigonometry  formula  gives  [7] 

sin  Az  =  -cosh  sin  h/cos<p.  (6) 
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Figure  2.  Relationship  between  xy  and  image  plane 
coordinate  systems. 


Figure  3  The  navigational  triangle 
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2.3  Image  point  position 

Let/  =  the  camera  lens  focal  length  =  distance  from 
I  lens  to  image  plane,  since  the  camera  is  assumed  to  be 

focused  at  infinity.  Then  it  is  easy  to  see  that 


allows  Equation  (12)  to  be  put  in  terms  of  a  function  of 
fa),  instead  of  2fa).  It  should  be  noted  that  only  cosu  is 
needed  for  the  subsequent  analysis,  so  it  is  not  necessary 
to  invert  the  above  equation  to  find  fa)  explicitly. 


tann  =  (x,2  +  y,2)  1/2 If,  (7) 

where  (x2,  yt)  are  the  convenient  set  of  image  plane  coor- 
>  dinates.  However,  v  is  referenced  to  the  x-y  axis  system, 

so  it  is  necessary  to  perform  the  transformation  between 
the  two  systems.  Reference  to  Figure  2  shows  that 


v  m  Az  +  Hdg  +  90°  -  tan-1  (yt/xt).  (8) 

2.4  Angle  of  incidence  or  reflection 

Consider  a  vector  A  that  points  from  the  reflecting  sur¬ 
face  toward  the  sun  (opposite  to  the  direction  of  the  inci¬ 
dent  ray)  and  another  vector  B  that  points  along  the 
reflected  ray  (from  the  sea  surface  to  the  image).  Then 
the  angle  between  A  and  the  normal  to  the  surface  is  the 
incidence  angle,  and  the  angle  between  B  and  the  surface 
normal  is  the  reflection  angle.  From  the  geometrical  op¬ 
tics  law  of  reflection  X,  B,  and  the  surface  normal  all 
lie  in  the  same  plane,  and  the  angle  of  incidence  equals 
the  angle  of  reflection  [5].  Let  this  angle  be  called  fa),  as 
defined  in  Section  2.1.  Then 

A  •  fT=  |iT|  |U|  cos  2fa),  (9) 

since2fa)  is  the  angle  between  X  and  B  If  the  components 
of  A  and  B  are  written  in  the  Xj-X2-X3  (sea  surface) 
coordinate  system, 

A,  =  0  Bt  =*  -|2f|  sinfi  sinO 

A  2  =  |  A  \  cos<f>  B2  =  |f?|  sinfi  cosd  (10) 

A5  =  |A|  sin<)>  B3  =  |2?|  cosu 

where  6  is  the  angle  between  the  sun  and  the  image  point 
P,  measured  counterclockwise  from  the  sun.  Since  the  im¬ 
age  coordinates  x  and  y  are  reflections  of  X1  and  X2, 
respectively,  it  is  easily  seen  that 

0  =  180°  -  v  (11) 

Then  the  result  of  performing  the  operations  implied  by 
Equation  (9)  is 


cos  2w  =  cos<(>  sinfi  cosd  +  sin<f>  cosfi 

*  -cos<(>  sinfi  cosy  +  sin<f>  cosfi. 

Application  of  the  trigonometric  identity 

cos2u  «  ‘si  (cos  2u  +  1) 


(12) 


2.5  Sea-surface  tilt 

The  sea-surface  tilt  angle  fi  can  be  found  from  Equa¬ 
tion  (2),  which  can  be  rewritten 

cosfi  =  (cosu  +  sin<f>)/2  cosu.  (14) 

2.6  Direction  of  steepest  ascent  of 
sea  surface 

From  the  foregoing  we  have  two  equations  for  the  local 
azimuth  of  ascent  of  sea  surface  a.  Equations  (1)  and  (3). 
Appendix  A  shows  how  Equation  (3)  may  be  solved.  Let 

A  =  1/tan  v,  (15a) 

B  =  cos<f)/2  sinfi  cosu.  (15b) 

Then 


a  -  tan~ 


AB±(A2+1-B2)» 
-B±A(A2+  1  -B2)‘A 


(16) 


where  the  ±  signs  are  paired,  +  with  +  or  -  with  - . 
There  are  two  solutions  for  a,  both  of  which  can  be  proved 
to  satisfy  Equation  (3).  Equation  (1)  gives  a  method  of  solv¬ 
ing  for  cosor,  again  there  are  two  solutions  for  a,  since 
cos  (-a)  —  cosa.  The  two  equations,  (1)  and  (3),  resolve 
the  ambiguity.  One  of  the  solutions  of  Equation  (3),  given 
by  Equation  (16),  matches  one  of  the  solutions  of  Equa¬ 
tion  (1);  that  one  is  the  correct  solution. 

It  should  be  noted  that  Cox  and  Munk  [2],  as  well  as 
Kinsman  [3],  show  Equation  (3)  with  a  minus  sign  affixed 
to  the  last  term  on  the  right  side.  This  is  not  borne  out 
by  the  derivation.  Furthermore,  if  the  equation  as  given 
in  those  references  is  solved,  giving  a  pair  of  solutions 
similar  to  those  given  by  Equation  (16),  it  is  found  that 
neither  of  those  solutions  satisfies  Equation  (1).  So,  clear¬ 
ly,  there  is  a  sign  error  in  those  references. 

Equation  (1)  relates  a  to  fi,  u,  and  <f>.  Equation  (3)  also 
relates  a  to  fi,  u,  and  <f>,  and  to  v  as  well.  This  may  be 
looked  at  in  two  ways: 

•  Equation  (3)  has  two  solutions  for  a  in  general,  only 
one  of  which  is  consistent  with  the  physics  of  the  prob¬ 
lem.  For  some  choices  of  fi,  u ,  0,  and  v.  Appendix 
A  shows  there  are  no  solutions  for  a.  But  in  this 
analysis  those  parameters  are  not  all  free.  The  angles 
fi,  fa),  and  <(>,  along  with  an  a  that  satisfies  Equation 
(3),  are  constrained  by  Equation  (1). 


(13) 


•  For  a  set  of  values  a,  0,  (*>,  and  <t>  that  satisfy  Equa¬ 
tion  (1),  Equation  (3)  fixes  the  value  of  v  (within  180°). 
But  v  is  given  by  Equation  (8). 

This  last  observation  makes  it  clear  that  all  of  the  angles 
are  functionally  related  in  a  complicated  way.  Equations 
(1),  (3),  (12),  and  (14)  (or  (2))  relate  a,  0 ,  /*,  v,  o>,  and 
<t>.  But  <f>,  n,  and  v  are  given  in  terms  of  external  condi¬ 
tions  by  Equations  (5),  (7),  and  (8),  respectively.  So  this 
leaves  three  variables  related  by  four  equations,  apparent¬ 
ly  one  more  than  necessary.  But  the  equations  all  involve 
trigonometric  functions  of  the  angles.  All  the  inverse 
trigonometric  functions  have  two  solutions  (apart  from 
the  periodicity  outside  the  range  -  180°  to  +  180°).  The 
physics  of  the  problem  resolves  much  of  the  ambiguity. 
The  angles  0,  /t,  to,  and  <f>  are  restricted  to  the  first 
quadrant.  The  value  of  v  can  be  found  uniquely  from  Equa¬ 
tion  (8);  since  both  x,  and  y [  are  known  separately,  the 
inverse  tangent  is  not  ambiguous.  But  a  may  take  on  any 
value  - 180°  <  a  <  - 180°,  and  it  is  not  immediately  clear 
how  to  eliminate  one  of  the  solutions  of  Equation  (1)  or 
(3)  without  using  the  other  equation. 

The  solution  process  for  Equation  (3)  derived  in  Ap¬ 
pendix  A  gives  an  expression  for  cos  a: 

-B±A(A2  +  1-B2)» 

cosa  =  - 1 - —  .  (17a) 

A2+  1 

This  may  be  equated  to  the  expression  for  cosa  that  is 
found  from  Equation  (1)  to  give  an  identity.  The  identity  is 

±(4  sin20  cos2(j)  -  sin2v  cos2<f)  )Vl 

=  sign  (sinv)  x  (sin/i  -  cosv  costj) ),  (17b) 

in  terms  of  the  original  variables.  The  ±  sign  on  the  left 
side  is  the  same  as  that  in  Equation  (16).  So  the  sign  of 
the  expression  on  the  right  side  of  Equation  (17b)  tells 
which  sign  to  pick  in  Equation  (16),  thus  resolving  the 
ambiguity. 

2.7  Distribution  over  a  and  0 

The  histogram  that  approximates  the  bivariate  distribu¬ 
tion  as  a  function  of  a  and  0  is  very  easily  found  by  set¬ 
ting  up  an  array,  each  of  whose  cells  (or  “bins”)  is  in¬ 
dexed  by  a  specific  range  of  a  and  0,  and  accumulating 
the  corresponding  pixel  values  in  each  cell.  The  choice 
of  those  “class  intervals”  for  a  and  0  is  of  some  con¬ 
cern.  Because  the  histogram  is  constructed  from  a  finite 
number  of  samples,  if  the  subdivision  is  too  fine  the  shape 
of  the  distribution  will  be  masked  by  statistical  fluctua¬ 
tions.  On  the  other  hand,  if  the  subdivision  is  too  coarse, 
the  shape  of  the  distribution  will  be  blurred  (an  extreme 


case  would  be  a  histogram  with  only  a  single  bin  that  con¬ 
tains  all  the  points).  Between  the  extremes  there  is  a  fair¬ 
ly  wide  range  of  reasonable  values. 

Sturges  [8]  has  given  a  criterion  for  choosing  a  class 
interval.  For  a  statistical  series  of  range  R  with  N  items, 
the  Sturges  criterion  for  the  optimal  class  interval  is 

C  =  R/(l  +  3.322  log  N).  (18) 

It  is  based  on  the  principle  that  the  proper  distribution 
into  classes  is  given,  for  all  numbers  that  are  powers  of 
2,  by  a  series  of  binomial  coefficients.  (Note  that  3-322 
=  l/log}02.)  For  example,  16  items  would  be  divided 
normally  into  5  classes,  with  class  frequencies  1,  4,  6,  4, 
1 .  So  if  a  statistical  series  had  16  items  with  values  rang¬ 
ing  from  20  to  70,  or  a  range  of  50  points,  it  should  be 
divided  into  five  classes  of  10  points  each;  the  class  inter¬ 
val  would  be  10. 

In  the  present  case,  each  digitized  image  consists  of  512 
x  512  =  262,144  =  2,H  pixels.  Thus,  there  are  262,144 
values  of  a  and  262,144  values  of  0.  From  Equation  (18), 
C  =  R/l  9.  Now,  in  practice,  the  class  interval  is  chosen 
to  be  some  convenient  value  near  the  theoretically  op¬ 
timum  interval.  If  we  use  C  =  R/l  8,  then  C  =  360° VI 8 
=  20°  for  a  and  C  =  90°/18  =  5°  for  0.  However, 
it  is  clear  that  0  =  90°  is  physically  impossible,  and  very 
large  (near  90°)  values  are  extremely  unlikely.  In  fact,  Cox 
and  Munk  report  that  0  =  25°  -  30°  may  be  a  prac¬ 
tical  maximum  [2],  So  the  range  R  within  which  values 
occur  is  less  than  90°;  consequently,  it  is  reasonable  to 
choose  a  class  interval  or  “bin  size”  smaller  than  5°.  An 
interval  of  1°  was  chosen  for  0.  For  a,  an  interval  of  10° 
(instead  of  20°)  was  chosen  to  obtain  better  resolution. 
“Bins”  of  these  widths  were  set  up  to  cover 
- 180°  <  a  ^  180°  and  0°  ^ 0  <90°.  So  the  histogram  ar¬ 
ray  contained  36  x  90  =  3240  cells. 

2.8  Distribution  over  components 
of  slope 

Cox  and  Munk  [1-3]  related  the  statistics  of  the  sea  sur¬ 
face  to  the  statistics  of  the  distribution  of  brightness  over 
components  of  slope.  Those  components  are 

Zx  =  tan0  sina*  (19a) 

Zy  =  tan0  cosa*  (19b) 

where  0  =  sea-surface  tilt  as  defined  before,  and  a  *  = 
azimuth  of  ascent  measured  clockwise  from  some  axis. 
When  that  axis  points  toward  the  sun,  a  *  =  a.  But  Cox 
and  Munk  found  that  the  principal  axes  of  the  distribu¬ 
tion  are  in  the  direction  of  the  wind  and  crosswind  [1-2]. 
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So  when  ct  *  is  measured  relative  to  the  wind  direction, 
the  principal  axes  are  the  coordinate  axes.  It  is  easily  shown 
that  this  a*  is  given  by 

a*  »  a  +  ( 180 0  -  Q  -  Ax)  (20) 

where  0  =  direction  of  wind,  measured  clockwise  from 
N  (similarly  to  Hdg)  and  a  and  Az  were  defined  previous¬ 
ly.  The  angle  Q  gives  the  direction  from  which  the  wind 
is  blowing,  in  agreement  with  common  usage:  a  west  wind 
is  one  that  blows  from  the  west. 

It  is  clear  that  [  - 1, 1]  is  likely  to  be  a  sufficient  range 
for  both  Zx  and  Z  ,  although  in  principle  they  are  un¬ 
bounded.  This  is  so  because  |.r/#Qf*j<I  and  jco.ra*|  ^1, 
and  even  when  |r/»a  *  |  =  1  or  | cosa  *  j  ==  1,  tanfi  ^  1  ex¬ 
cept  for  /3  >  45  °,  which  is  unlikely.  So  the  histogram  was 
set  up  for  - 1  ^  Zx  ^  1  and  - 1  ^  Zy  ^  1.  The  a-|3 
histogram  has  36  x  45  =  1620  cells  over  the  range  of 
(8  up  to  45°.  The  Zx-Zy  histogram  was  given  approx¬ 
imately  the  same  number  of  cells  by  dividing  the  Zx  and 
Zy  ranges  into  41  intervals  each  (the  odd  number  causes 
one  interval— the  central  one— to  be  centered  at  zero). 
The  histogram  is  constructed  by  calculating  Zx  and  Zy 
for  each  pixel,  then  adding  the  pixel  value  to  the  ap¬ 
propriate  cell. 

2.9  Correction  for  nonzero  roll 

The  equations  for  the  sun  glitter  analysis  that  were  de¬ 
rived  in  preceding  sections  are  based  on  a  geometry  in 
which  the  photographic  film  is  horizontal.  Figure  4 


Figure  4.  Effect  of  camera  tilt  on  image  point  position. 


illustrates  that  when  the  film  is  not  horizontal  the  intercept 
of  a  light  ray  with  the  film  changes,  in  general.  The  axis 
of  rotation,  which  passes  through  the  center  of  the 
camera’s  aperture,  is  normal  to  the  plane  of  the  figure. 
The  dashed  lines  represent  the  axis  of  the  camera’s  imag¬ 
ing  geometry,  normal  to  the  film  in  each  case.  The  solid 
arrow  shows  a  light  ray,  which  intersects  the  horizontal 
and  tilted  image  planes  in  different  positions  in  the  respec¬ 
tive  coordinate  systems.  It  is  clear  that  there  is  a  similar 
projective  effect  on  the  component  perpendicular  to  the 
plane  of  the  figure  (except  along  the  line  of  intersection 
of  the  two  film  planes). 

The  effect  enters  the  analysis  through  p  and  i>,  and 
alters  the  values  found  for  w,  j3,  and  a.  If  the  plane  of 
the  film  is  tilted  and  p  and  v  are  calculated  from  the  (tilted) 
image  point  position  by  the  procedure  of  Section  2.3,  the 
a  and  /3  that  result  will  be  in  error. 

A  solution  is  to  calculate  a  corrected  set  of  p,  v  that 
correspond  to  the  position  on  a  hypothetical  horizontal 
film  that  the  same  light  ray  would  reach.  To  do  this,  define 
two  coordinate  systems  xJy1-zi  and  x2y2-z2  with  a  com¬ 
mon  origin  O  at  the  camera's  aperture.  The  former  system 
has  the  xly1  plane  horizontal  with  zt  pointing  upward. 
The  x2-y2-z2  system  is  rotated  an  angle  r  about  the  com¬ 
mon  y  axis  with  respect  to  the  other  coordinate  system; 
r  is  positive  for  a  clockwise  rotation  from  system  1  to 
system  2. 

The  x2-y2  z2  system  is  oriented  according  to  the  at 
titude  of  the  aircraft  with  nonzero  roll  angle  r  the  ^-axis 
points  toward  the  starboard  wing  and  the  jv^-axis  points 
toward  the  nose.  Figure  5  illustrates  the  two  coordinate 
systems. 

The  (hypothetical)  film  plane  in  the  first  system  is  the 
plane  z  j  =  /;  in  the  second  system  the  (true)  film  plane 
is  z2  =  f.  Here  /is  the  camera’s  focal  length.  The  xy 
coordinates  of  points  in  the  film  plane  (either  one)  will 
be  denoted  by  an  extra  subscript  “f.”  Thus,  (x,^  yL)  cor¬ 
respond  to  the  ( x j,  yt)  that  were  used  in  Section  2.3  to 
calculate  p  and  v.  They  are  the  quantities  needed  here. 

The  transformation  between  arbitrary  coordinates  in 


Now  consider  a  light  ray  that  passes  through  O  and  strikes 
the  actual  film  plane  z2  =  /  at  the  point  P2  -  (z2  y^ 
f).  We  wish  to  find  the  point  where  the  same  ray  would 
strike  the  hypothetical  film  plane  z1  =  f.  This  can  be  easi¬ 
ly  done  in  three  stages.  First,  from  Equation  (21),  the 
coordinates  of  P2  in  the  xty,z,  reference  frame  are 


5 


From  this  we  can  calculate  values  for  /x  and  v  consistent 
with  the  geometry  that  was  originally  assumed. 

From  Equation  (7) 

tan  it.  =  (x,f  +  y,/)'A/f 

(x2/+y2/)+fsin  r [f  sin  r+x2/(cos  r  +  D)]\A 


where  u  -  cos  r  -  x^sin  r/f.  Squaring  Equation  (25), 


Figure  5.  Relationship  between  horizontal  and  rolled  coor¬ 
dinate  system. 


x2fcos  r  +  /  sin  r 


yi\  = 


-x 2f  sin  r  +  f  cos  r 


These  are  not  the  coordinates  of  the  point  Pt  where  the 
light  ray  strikes  the  hypothetical  film  plane  z1  =  /.To 
find  those  coordinates  we  next  write  the  equation  of  the  line 
that  passes  through  the  points  O  and  P2.  Since  O  is  the 
origin,  substituting  from  Equation  (22),  that  equation  is 


x 2f  cos  r  +  f  sin  r  y2j  f  cos  r  -  x2j  sin  r 

(23) 

This  is  the  equation  of  the  light  ray  in  the  x1-yl-z1 
system.  From  it,  the  coordinates  of  P,  are 


/  (f  sin  r  +  x2j  cos  r) 
{f  cos  r  -  x2fsin  r) 


y,f  = 


( f  cos  r  -  x2f  sin  r) 


Z  is  —  f 


sin  r  if  sin  r+x?f(cos  r  +  D)l  I 

- - -2- - i  j*  .  (26) 

The  value  of  v  is  found  from  Equation  (8),  where  Equa¬ 
tions  (24a)  and  (24b)  must  be  substituted  for  x,  and  yh 
respectively. 

3.  Limitations  of  analysis 

The  foregoing  text  has  described  the  derivation  of  the 
sun  glitter  analysis  algorithms.  The  discussion  will  con¬ 
clude  with  a  summary  of  the  approximations  involved, 
some  of  which  are  not  explicit. 

Three  of  the  approximations  relate  to  the  geometry  of 
the  problem: 

•  A  spherical  earth  was  assumed. 

•  No  refraction  correction  was  applied  to  the  apparent 
position  of  the  sun. 

•  Mean  solar  time  was  used  in  calculating  the  sun’s 
position. 

All  have  some  effect  on  the  sun’s  position  in  the  assumed 
geometry.  The  first  causes  relatively  negligible  effects.  The 
second  causes  an  error  of  up  to  about  Vi0  in  the  sun's 
elevation.  The  error  due  to  the  third  can  cause  an  error 
of  up  to  4°  in  the  local  hour  angle.  Both  errors  are  then 
propagated  to  the  quantities  that  depend  on  <)>  and  h. 
Two  other  approximations  involve  photometry: 

•  No  correction  of  the  digitized  values  to  obtain  original 
luminance  values  was  performed. 

•  No  background  correction  was  performed. 

With  respect  to  the  first  of  these,  the  digitized  positive 
image  brightness  values  are  almost  certainly  not  directly 
proportional  to  the  reflected  brightness  field  that  il¬ 
luminated  the  film.  It  is  the  latter  that  should  be  used 


to  construct  the  brightness  distributions.  While  the  digit¬ 
ized  values  are  probably  related  to  the  original  brightness 
values  by  a  monotonically  increasing  function,  the  form 
of  that  function  is  unknown.  It  involves  the  properties 
of  the  camera's  lens,  the  response  of  the  film  (D  -  log 
E  curve),  the  uniformity  of  illumination  of  the  transparency 
when  it  is  digitized,  the  response  of  the  digitizing  camera, 
and  the  law  by  which  that  camera's  output  is  transformed 
to  digitized  values.  The  background  correction  refers  to 
compensation  for  the  radiation  from  the  reflection  of 
skylight  at  the  sea  surface  and  from  sunlight  scattered  by 
particles  beneath  the  sea  surface.  It  was  implicitly  assumed 
in  the  foregoing  discussion  that  all  of  the  light  that  reaches 
the  film  is  from  specular  reflection  of  sunlight  from  the 


sea  surface.  Cox  and  Munk  devoted  considerable  atten 
tion  to  background  correction  [2]. 


4.  Program  documentation 

This  section  contains  formal  documentation  of  a  main 
program  and  15  subroutines  that  were  written  (in  FOR 
TRAN  IV)  to  perform  sun  glitter  analysis  of  aerial 
photographs.  Each  routine's  documentation  is  self- 
contained.  Thus,  in  a  few  cases,  equations  are  numbered 
independently  starting  with  Equation  (1).  The  written 
documentation  is  followed  by  user  instructions,  program 
listings,  and  a  sample  run  stream. 
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1.  Name 


GMAIN2 

2.  Purpose 

The  purpose  of  GMAIN2  is  to  calculate  histograms 
that  approximate  distributions  of  sea-surface  slope  by 
analyzing  aerial  photographs  of  the  sun’s  reflection  on 
the  water. 

3.  Calling  sequence 
GMAIN2  is  a  main  program. 

4.  Input — output 

4.1  Input 

The  following  items  are  read  from  logical  unit  5  (user 


terminal): 

LAT 

latitude 

(real) 

LONG 

= 

longitude 

(real) 

DECL 

= 

sun’s  declination 

(real) 

TIME 

= 

Zulu  time 

(real) 

HDG 

= 

aircraft  heading 

(real) 

ROLL 

= 

aircraft  roll  angle 

(real) 

WIND 

= 

wind  direction 

(real) 

F 

= 

camera  focal  length 

(real) 

HEIGHT 

= 

image  height  (or  width) 

(real) 

IR1 

= 

first  row  of  region  to  process 

(integer) 

IR2 

= 

last  row  of  region  to  process 

(integer) 

IE1 

= 

first  column  of  region  to  process 

(integer) 

IE2 

= 

last  column  of  region  to  process 

(integer) 

FILE 

= 

image  file  name  (CHARACTER'S) 

The  above  items  are  read  once  per  run.  All  angles  are 
in  degrees. 

In  addition,  records  containing  image  data  are  read  from 
logical  unit  1.  Each  record  consists  of 

BUFFER  =  312-word  array  that  contains  two  con¬ 
secutive  scan  lines,  1  pixel  per  byte.  (integer) 

4.2  Output 

Messages  that  constitute  the  other  half  of  an  interac¬ 
tive  dialogue  are  displayed  on  logical  unit  6  (user  terminal). 
The  messages  solicit  the  items  that  are  read  from  logical 
unit  5,  as  described  in  4.2. 

All  items  read  from  logical  unit  3  are  printed  on  logical 
unit  7  (line  printer).  In  addition,  the  following  items  are 
printed. 


AZ  =  solar  azimuth  (degrees)  (real) 

HIST  =  distribution  of  integrated  image  (real) 

intensity  vs.  a  and  0 

A1  =  minimum  a  (real) 

A2  =  maximum  a  (real) 

DALPHA  =  a  increment  (real) 

B1  =  minimum  0  (real) 

B2  =  maximum  0  (real) 

DBETA  =  0  increment  (real) 

HISTB  =  HIST  summed  over  a,  as  a  (real) 

function  of  0 

XMEANB  =  mean  of  0  with  respect  to  (real) 


distribution  HISTB 

VARB  =  variance  of  0  with  respect  to  (real) 

distribution  HISTB 

HISTCU  =  distribution  of  integrated  image  (real) 

intensity  vs.  crosswind  and  up¬ 
wind  components  of  sea-surface 
slope 

ZMIN  =  minimum  value  of  slope  components  (real) 

ZMAX  =  maximum  value  of  slope  components  (real) 

ZDELT  =  slope  component  increment  (real) 

SLICE  =  first  crosswind,  then  upwind  slice  (real) 

(along  coordinate  axes)  of 
HISTCU(both  are  printed) 

RMS  =  root-mean-square  value  of  slope  (real) 

in  slice  (computed  for  both) 

In  the  event  of  certain  error  conditions  in  the  calculations 
performed  by  subroutines,  appropriate  messages  are  printed 
on  logical  unit  7. 

4.3  File  Storage 

None. 

3.  Exits 

There  are  no  nonstandard  exits. 

6.  Usage 

GMAIN2  is  written  in  FORTRAN  IV  and  is  present¬ 
ly  implemented  on  the  HP- 3000  operating  under  MPE-III. 

7.  External  interfaces 

7.1  System  subroutines 

DCOS 

DSIN 

7.2  Other  programs  called 

GETPIX 

HDSPLY 


MNVAR 

NRMLZ 

RMSARY 

STEP1 

STEP2 

STEP3 

STEP4 

STEP5 

STEP6 

STEP? 

7.3  External  storage  used 
None. 

8.  Performance  specifications 

8.1  Storage 

Stack:  11,765  words 
Code:  2686  words 

8.2  Execution  time 
TBD. 

8.3  I/O  Load 

I/O  statements  are  described  in  4.1  and  4.2.  The  number 
of  records  read  from  unit  1  depends  on  IR1,  IR2,  IE  1,  IE2. 


8.4  Restrictions 

Image  files  are  assumed  to  be  in  the  format  produced 
by  Cl. 

9.  Method 

Constants  are  read  in  and  program  variables  are  initializ¬ 
ed.  GMAIN2  then  skips  to  the  first  record  needed  from 
logical  unit  1.  STEP1  is  called  to  calculate  solar  elevation 
and  azimuth,  and  the  correction  to  reference  a  to  the 
wind  direction  is  computed.  The  first  record  is  read  from 
unit  1  and  unpacked  to  one  pixel  per  word.  Then  a  loop 
begins  over  the  rows  and  columns  to  be  processed.  Coor¬ 
dinates  of  the  current  image  point  are  calculated,  and 
STEP2-STEP5  are  called  to  calculate  (ultimately)  the  a 
and  0  for  that  point.  STEP6  and  STEP7  are  then  called 
to  increment  HIST  and  HISTCU.  At  the  end  of  the  loop, 
another  record  is  read  from  unit  1  when  necessary.  The 
remainder  of  the  program  is  concerned  with  the  calcula¬ 
tion  of  statistical  measures  and  with  output. 

10.  Comments 

GMAIN2  does  not  include  a  calculation  to  recover  the 
original  light  intensity  from  digitized  image  values. 


1.  Name 


7.  External  interfaces 


STEP1 

2.  Purpose 

The  purpose  of  STEP1  is  to  calculate  local  hour  angle, 
solar  elevation,  and  solar  azimuth.  The  last  two  are  return¬ 
ed  to  the  calling  routine. 

3.  Calling  sequence 

CALL  STEP1  (TIME,  LONG,  LAT,  DECL,  PHI, 


CDPHI,  SDPHI,  AZ,  $n). 
where: 

TIME  =  Zulu  time  (real) 

LONG  =  longitude  W  of  Greenwich  (real) 

LAT  =  latitude  N  of  equator  (real) 

DECL  =  solar  declination  (real) 

PHI  =  solar  elevation  (real) 

CDPHI  =  cosine  of  PHI  (double  precision) 

SDPHI  =  sine  of  PHI  (double  precision) 

AZ  =  solar  azimuth  (real) 


n  =  statement  number  in  calling  program 

to  which  control  is  returned  if  AZ 
cannot  be  calculated 

4.  Input— output 

4.1  Input 

There  are  no  input  statements. 

4.2  Output 

There  are  no  output  statements. 

4.3  File  Storage 
None. 

5.  Exits 

If  the  expression  for  (sin  (AZ)|  >  1,  AZ  is  set  to  zero 
and  the  error  return  is  taken. 

6.  Usage 

STEP1  is  written  in  FORTRAN  IV  and  is  presently 
implemented  on  the  HP-3000  operating  under  MPE-III. 


7.1  System  subroutines 

DCOS 

DSIN 

7.2  Other  programs  called 
ASIN 

7.3  External  storage  used 
None. 

8.  Performance  specifications 

8.1  Storage 

Stack:  56  words 
Code:  231  words 

8.2  Execution  time 
TBD. 

8.3  I/O  Load 

There  are  no  I/O  statements. 

8.4  Restrictions 
None. 

9.  Method 

Hour  angle  is  calculated  by 
h  =  {(t  -  720)/4j  -  X  . 

Solar  elevation  is  calculated  by 
sin<t>  =  sinb  sin  L  +  cosb  cos  L  cos  h 
Solar  azimuth  is  calculated  by 
sin  Az  =  -cosb  sin  h/cos(f>  . 

All  symbols  are  defined  in  Section  2  of  this  report. 

10.  Comments 
None. 
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1.  Name 


7.  External  interfaces 


STEP2 

Source  file  STEP2A 

2.  Purpose 

The  purpose  of  STEP2  is  to  calculate  image  angular 
coordinates. 

3.  Calling  sequence 

CALL  STEP2  (X,  Y,  F,  HDG,  AZ,  MU,  NU,  CDR, 
SDR,  IFLAG). 

where: 

X  =  image  coordinate  toward  starboard  wing  (real) 

Y  =  image  coordinate  toward  aircraft  nose  (real) 


F 

=  camera  lens  focal  length 

(real) 

HDG 

*=  aircraft  heading 

(real) 

AZ 

=  solar  azimuth 

(real) 

MU 

=  angle  between  reflected  ray  and  vertical  (real) 

NU 

=  angle  of  image  point  measured 

(real) 

clockwise  from  an  axis  that  points 
away  from  the  sun 

CDR  =  cosine  of  aircraft  roll  angle  (double  precision) 
SDR  =  sine  of  aircraft  roll  angle  (double  precision) 
IFLAG  =  flag  that  indicates  error  in  MU  (integer) 
calculation  when  roll  is  nonzero 
=  0— no  error 
=■>  1 — error 

4.  Input — output 

4.1  Input 

There  are  no  input  statements. 

4.2  Output 

There  are  no  output  statements. 

4.3  File  Storage 
None. 

5.  Exits 

There  are  no  nonstandard  exits. 

6.  Usage 

STEP2  is  written  in  FORTRAN  IV  and  is  presently 
implemented  on  the  HP-3000  operating  under  MPE-DI. 


7.1  System  subroutines 

DATAN 

DATAN2 

DSQRT 

7.2  Other  programs  called 
None. 

7.3  External  storage  used 
None. 

8.  Performance  specifications 

8.1  Storage 

Stack:  43  words 
Code:  320  words 

8.2  Execution  time 
TBD. 

8.3  I/O  Load 

There  are  no  I/O  statements. 

8.4  Restrictions 
None. 

9.  Method 

If  the  roll  angle  r  =  0,  the  output  quantities  are 
calculated  from 

tanfi  =  (xj2  +  yt2),/2/f>  and 

■  Az  +  Hdg  +  90°  -tan~J  (yt/xt) 

where  the  symbols  are  defined  in  Section  2  of  this  report. 
In  the  above  equations  Xj  =  X  and  yt  =  Y,  cartesian 
coordinates  in  the  image. 

If  r  ^  0,  then 

tann  =  (xlf*  +  ytf)*/f 

_  {  (*#*+y#*)+fri"  '[/***  r+x2f(cosr  +  EHil 


1.  Name 


7.  External  interfaces 


STEP3 

2.  Purpose 

The  purpose  of  STEP3  is  to  calculate  the  cosine  of  the 
angle  of  incidence  or  reflection. 

3.  Calling  sequence 

CALL  STEP3  (CDPHI,  SDPHI,  MU,  NU,  COSOMG). 
where: 

CDPHI  =  cosine  of  solar  elevation  (double  precision) 

SDPHI  =  sine  of  solar  elevation  (double  precision) 

MU  **  angle  between  reflected  (real) 

ray  and  vertical 

NU  =  angle  of  image  point  measured  (real) 

clockwise  from  an  angle  that 
points  away  from  the  sun 

COSOMG  =  cosine  of  angle  of  (double  precision) 
incidence  or  reflection 

4.  Input— output 

4.1  Input 

There  are  no  input  statements. 

4.2  Output 

There  are  no  output  statements. 

4.3  File  Storage 
None. 

5.  Exits 

There  are  no  nonstandard  exits. 

6.  Usage 

STEP3  is  written  in  FORTRAN  IV  and  is  presently 
implemented  on  the  HP-3000  operating  under  MPE-III. 


7.1  System  subroutines 

DCOS 

DSIN 

DSQRT 

7.2  Other  programs  called 
None. 

7.3  External  storage  used 
None. 

8.  Performance  specifications 

8.1  Storage 

Stack:  26  words 
Code:  133  words 

8.2  Execution  time 
TBD. 

8.3  I/O  Load 

There  are  no  I/O  statements. 

8.4  Restrictions 
None. 

9.  Method 

The  equations  that  are  implemented  are 

cos  2o}  =  -cos<l>  sinfi  cosv  +  sitt<f>  cosfi, 

cosos  =  {‘A  (cos  *  jj  +  1)}H- 

The  symbols  are  defined  in  Section  2  of  this  report.  It 
is  assumed  that  cu  is  in  the  first  quadrant. 

10.  Comments 
None. 


1.  Name 


7.  External  interfaces 


STEP4 

2.  Purpose 

The  purpose  of  STEP4  is  to  calculate  sea-surface  tilt. 

3.  Calling  sequence 

CALL  STEP4  (MU,  SDPHI,  COSOMG,  BETA) 
where: 

MU  =  angle  between  reflected  ray  (real) 

and  vertical 

SDPHI  =  sine  of  solar  elevation  (double  precision) 

COSOMG  =  cosine  of  angle  of  (double  precision) 

incidence  or  reflection 

BETA  =  sea-surface  tilt  (real) 

4.  Input— output 

4.1  Input 

There  are  no  input  statements. 

4.2  Output 

There  are  no  output  statements. 

4.3  File  Storage 
None. 

3.  Exits 

There  are  no  nonstandard  exits. 

6.  Usage 

STEP4  is  written  in  FORTRAN  IV  and  is  presently 
implemented  on  the  HP- 3000  operating  under  MPE-III. 


7.1  System  subroutines 
DCOS 

7.2  Other  programs  called 
ACOS 

7.3  External  storage  used 
None. 

8.  Performance  specifications 

8.1  Storage 

Stack:  14  words 
Code:  64  words 

8.2  Execution  time 
TBD. 

8.3  I/O  Load 

There  are  no  I/O  statements. 

8.4  Restrictions 
None. 

9-  Method 

The  sea-surface  tilt  angle  is  found  from 

cos0  =  (cosh  +  sin<f>)/2  cosoj, 

where  the  symbols  are  defined  in  Section  2  of  this  report. 
If  cos  «  =  0,  which  implies  glancing  incidence,  /3  is  set 
to  zero. 

10.  Comments 

None. 


1.  Name 


root  is  set  to  zero  in  the  ALPHA  calculation  and  the  er¬ 
ror  return  is  taken. 


STEP5 

Source  file  STEP5A 

2.  Purpose 

The  purpose  of  STEP5  is  to  calculate  the  local  azimuth 
of  ascent  of  the  sea  surface. 

3.  Calling  sequence 

CALL  STEP5  (NU,  BETA,  PHI,  CDPHI,  COSOMG, 


MU,  ALPHA,  $n). 

where: 

NU 

= 

angle  of  image  point  measured 
clockwise 

(real) 

BETA 

= 

sea-surface  tilt 

(real) 

PHI 

= 

solar  elevation 

(real) 

CDPHI 

= 

cosine  of  PHI  (double  precision) 

COSOMG 

= 

cosine  of  angle  of  (double  precision) 

incidence  or  reflection 

MU 

— 

angle  between  reflected  ray  and 
vertical 

(real) 

ALPHA 

* 

local  azimuth  of  ascent  of  sea 

(real) 

surface 

n  =  statement  number  in  calling 

program  to  which  control  is 
returned  if  ALPHA  cannot  be 
calculated 

4.  Input— output 

4.1  Input 

There  are  no  input  statements. 

4.2  Output 

There  are  no  output  statements. 

4.3  File  Storage 

None. 

5.  Exits 


6.  Usage 

STEP5  is  written  in  FORTRAN  IV  and  is  presently 
implemented  on  the  HP-3000  operating  under  MPE-IIL 

7.  External  interfaces 

7.1  System  subroutines 

DATAN2 

DCOS 

DSIN 

DSQRT 

DTAN 

7.2  Other  programs  called 
ACOS 

7.3  External  storage  used 
None. 

8.  Performance  specifications 

8.1  Storage 

Stack:  71  words 
Code:  359  words 

8.2  Execution  time 
TBD. 

8.3  I/O  Load 

There  are  no  I/O  statements. 

8.4  Restrictions 
None. 

9.  Method 

The  local  azimuth  of  sea-surface  ascent  is  found  from 

» .  r  tteuit’-rji-i 


In  case  the  square  root  in  the  expression  for  ALPHA 
(see  Section  9)  involves  a  negative  number,  that  square 


where 


A  =  1/tan  v, 

B  =  cos<f>/2  sin&  cosu ), 

and  the  other  symbols  are  defined  in  Section  2  of  this 
report,  which  also  shows  that  the  sign  to  use  in  both  the 
numerator  and  the  denominator  of  tana  is  the  sign  of 
the  expression 


Special  cases  are  handled  as  described  in  Appendix  A. 
When  A2  +  1  -  B2  <  0,  the  square  root  cannot  be 
calculated.  This  should  occur  only  as  a  result  of  roundoff 
error  when  A2  +  1  -  B2  is  very  near  zero,  so  in  this 
case  a  is  calculated  from 

a  =  tan~1(AB/-B) 

and  the  error  return  is  taken. 

10.  Comments 


sinv  (sinn  -  cosv  cos<f>). 


None. 


1.  Name 


6.  Usage 


STEP6 

2.  Purpose 

The  purpose  of  STEP6  is  to  increment  the  a-0 
histogram. 

3.  Calling  sequence 

CALL  STEP6  (ALPHA,  BETA,  POINT,  HIST, 
NALPHA,  NBETA,  DALPHA,  DBETA). 

where: 

ALPHA  =  azimuth  of  ascent  of  sea  surface  (real) 


BETA  =  sea-surface  tilt  (real) 

POINT  =  image  point  intensity  value  (integer) 

HIST  =  histogram  array  (real) 

NALPHA  =  number  of  ALPHA  increments  (integer) 

NBETA  =  number  of  BETA  increments  (integer) 

DALPHA  =  size  of  an  ALPHA  increment  (integer) 

DBETA  =  size  of  a  BETA  increment  (integer) 


4.  Input— output 

4.1  Input 

There  are  no  input  statements. 

4.2  Output 

There  are  no  output  statements. 

4.3  File  Storage 
None. 

5.  Exits 

There  are  no  nonstandard  exits. 


STEP6  is  written  in  FORTRAN  IV  and  is  presently 
implemented  on  the  HP-3000  operating  under  MPE-III. 

7.  External  interfaces 

7.1  System  subroutines 
None. 

7.2  Other  programs  called 
None. 

7.3  External  storage  used 
None. 

8.  Performance  specifications 

8.1  Storage 

Stack:  4  words 
Code:  48  words 

8.2  Execution  time 
TBD. 

8.3  I/O  Load 

There  are  no  I/O  statements. 

8.4  Restrictions 
None. 

9.  Method 

The  indices  of  the  appropriate  histogram  cell  are 
calculated  from  ALPHA  and  BETA,  and  the  pixel  value 
is  added  to  that  cell. 

10.  Comments 
None. 
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1.  Name 


7.  External  interfaces 


STEP7 

2.  Purpose 

The  purpose  of  STEP7  is  to  increment  the  crosswind- 
upwind  histogram. 

3.  Calling  sequence 

CALL  STEP7  (ALPHA,  BETA,  DELTA,  POINT, 
HIST,  N,  ZDELT,  ZMAX). 

where: 


ALPHA  =  azimuth  of  ascent  of  sea  surface  (real) 
BETA  =  sea-surface  tilt  (real) 

DELTA  =  correction  to  reference  ALPHA  (real) 
to  wind  direction 

POINT  =  image  point  intensity  value  (integer) 
HIST  =  histogram  array  (real) 

N  =  number  of  increments  along  (integer) 

each  histogram  axis 

ZDELT  =  size  of  an  increment  (real) 

ZMAX  =  upper  limit  for  each  axis  (real) 


4.  Input— output 

4.1  Input 

There  are  no  input  statements. 

4.2  Output 

There  are  no  output  statements. 

4.3  File  Storage 
None. 

3.  Exits 

There  are  no  nonstandard  exits. 

6.  Usage 

STEP7  is  written  in  FORTRAN  IV  and  is  presently 
implemented  on  the  HP-3000  operating  under  MPE-III, 


7.1  System  subroutines 

COS 

SIN 

TAN 

7.2  Other  programs  called 
None. 

7.3  External  storage  used 
None. 

8.  Performance  specifications 

8.1  Storage 

Stack:  14  words 
Code:  91  words 

8.2  Execution  time 
TBD. 

8.3  I/O  Load 

There  are  no  I/O  statements. 

8.4  Restrictions 
None. 

9.  Method 

DELTA  is  used  to  reference  ALPHA  to  the  wind  direc¬ 
tion.  Then  the  two  histogram  coordinates  are  found  from 

Zx  =  tanfisina", 

Z  =  tan&  cos  a  *, 

The  symbols  are  defined  in  Section  2  of  this  report.  Then 
the  indices  of  the  histogram  cell  that  contains  Zx  and  Z 
are  calculated,  and  the  pixel  value  is  added  to  that  cell. 
Since  tan/3  is  unbounded  in  principle,  the  cell  indices  are 
corrected  if  necessary  to  lie  within  the  assumed  range. 

10.  Comments 


None. 


1.  Name 


7.  External  interfaces 


GETPIX 
(GET  PIXels) 

2.  Purpose 

The  purpose  of  GETPIX  is  to  unpack  a  read  buffer  that 
contains  two  image  lines  packed  one  pixel  per  byte  into 
a  two-row  array  in  which  each  word  contains  a  pixel. 

3.  Calling  sequence 

CALL  GETPIX  (BUFFER,  ROW), 
where: 

BUFFER  =  read  buffer 
ROW  =  two-row  array 

4.  Input— output 

4.1  Input 

There  are  no  input  statements. 

4.2  Output 

There  are  no  output  statements. 

4.3  File  Storage 
None. 

5.  Exits 

There  are  no  nonstandard  exits. 

6.  Usage 

GETPIX  is  written  in  FORTRAN  IV  and  is  presently 
implemented  on  the  HP-3000  operating  under  MPE-III. 


7.1  System  subroutines 
None. 

7.2  Other  programs  called 
None. 

7.3  External  storage  used 
None. 

8.  Performance  specifications 
8.1  Storage 


8.2  Execution  time 
TBD. 

8.3  I/O  Load 

There  are  no  I/O  statements. 

8.4  Restrictions 
None. 

9.  Method 

The  first  256  words  of  BUFFER  are  unpacked  into  the 
first  row  of  ROW,  with  each  consecutive  byte  from  the 
former  going  in  order  into  words  in  the  latter.  The  last 
256  words  of  BUFFER  are  unpacked  similarly  into  the 
second  row  of  ROW. 

10.  Comments 
None. 


(integer)  Stack:  3  words 

(integer)  Code:  48  words 
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1.  Name 
HDSPLY 

2.  Purpose 

The  purpose  of  HDSPLY  is  to  display  a  univariate 
histogram  on  the  line  printer. 

3.  Calling  sequence 

CALL  HDSPLY  (HIST,  N,  NAME,  XO,  STEP). 

where: 

HIST  =  histogram  array  (real) 

N  =  number  of  histogram  cells  (integer) 

NAME  =  name  of  independent  (CHARACTER*  10) 

variable 

XO  =  lower  end  of  independent  variable  (real) 

range 

STEP  =  independent  variable  step  size  (real) 

4.  Input— output 

4.1  Input 

There  are  no  input  statements. 

4.2  Output 

A  heading  and  the  histogram  itself  are  printed  on  logical 
unit  7.  The  heading  is 

(value  of  NAME)  OCCUPANCY 

The  histogram  is  printed  by  giving  “breakpoint”  values 
of  the  independent  variable  and  cell  occupancies.  The  items 
printed  are: 

X  =  current  breakpoint  value  (real) 

HIST(I)  =  histogram  occupancy  in  cell  I  (real) 

STAR  =  ASCII  .  (CHARACTER'D 

X  and  HIST(I)  are  interleaved  in  such  a  way  that  each 
HIST(I)  is  printed  between  values  of  the  independent 
variable  that  mark  the  boundaries  of  that  cell.  Along  with 
each  HIST(I),  STAR  is  printed  a  number  of  times  that 
is  proportional  to  HIST(I). 

4.3  File  Storage 
None. 

3.  Exits 


6.  Usage 

HDSPLY  is  written  in  FORTRAN  IV  and  is  presently 
implemented  on  the  HP-3000  operating  under  MPE-III. 

7.  External  interfaces 

7.1  System  subroutines 
None. 

7.2  Other  programs  called 
None. 

7.3  External  storage  used 
None. 

8.  Performance  specifications 

8.1  Storage 

Stack:  11  words 
Code:  174  words 

8.2  Execution  time 
TBD. 

8.3  I/O  Load 

Each  call  of  HDSPLY  produces  2N  +  2  lines  of  out¬ 
put  on  logical  unit  7. 

8.4  Restrictions 
None. 

9.  Method 

HIST  is  searched  to  find  the  maximum  occupancy, 
which  is  used  to  calculate  a  scale  factor  such  that  a  row 
of  100  asterisks  is  printed  for  the  maximum-occupancy 
cell  and  a  proportionate  number  of  asterisks  for  each  other 
cell.  Next,  the  title  is  printed.  This  is  followed  by  lines 
of  independent  variable  values  alternating  with  lines  giv¬ 
ing  cell  occupancy  and  a  row  of  asterisks  for  illustration. 
Between  0  and  100  asterisks  are  printed.  See  4.2  for  more 
information  on  the  output. 

10.  Comments 


There  are  no  nonstandard  exits. 


None. 


1.  Name 


ACOS 

Source  file  MLACOS 
2.  Purpose 

The  purpose  of  ACOS  is  to  calculate  the  single-precision 
inverse  cosine  of  a  value. 

3-  Calling  sequence 

Y  =  ACOS  (X) 

where: 

ACOS  =  function  value  (radians)  (real) 

X  =  cosine  of  angle  to  be  determined  (real) 

4.  Input— output 

4.1  Input 

There  are  no  input  statements. 

4.2  Output 

There  are  no  output  statements. 

4.3  File  Storage 
None. 

Y  Exits 

There  are  no  nonstandard  exits. 

6.  Usage 

ACOS  is  written  in  FORTRAN  IV  and  is  presently 
implemented  on  the  HP-3000  operating  under  MPE-111. 

7.  External  interfaces 

7.1  System  subroutines 

None. 


7.2  Other  programs  called 
AS1N 

7.3  External  storage  used 
None. 

8.  Performance  specifications 

8.1  Storage 

Stack:  0  words 
Code:  11  words 

8.2  Execution  time 
TBD. 

8.3  I/O  Load 

There  are  no  I/O  statements. 

8.4  Restrictions 
None. 

9.  Method 

The  arccosine  is  calculated  from 

cos~l  x  =  7T/2  -  sin~Ix. 

This  method  is  used  in  Reference  9.  ACOS  returns  the 
principal  value  of  the  arccosine.  That  is,  if  y  =  cas~  lx. 

0  ^  1 r/2  for  x  ^  0 

7r/2  <y  ^  7T  for  x  <  0 

10.  Comments 

See  the  documentation  for  ASIN  in  order  to  complete 
the  above  discussion. 
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1.  Name 
ASIN 

Source  file  MLASIN 

2.  Purpose 

The  purpose  of  ASIN  is  to  calculate  the  single-precision 
inverse  sine  of  a  value. 

3.  Calling  sequence 
Y  =  ASIN  (X) 


where: 

ASIN  =  function  value  (radians) 

X  =  sine  of  angle  to  be  determined 

4.  Input— output 

4.1  Input 

There  are  no  input  statements. 

4.2  Output 

There  are  no  output  statements. 

4.3  File  Storage 


3.  Exits 

There  are  no  nonstandard  exits. 

6.  Usage 

ASIN  is  written  in  FORTRAN  IV  and  is  presently  im¬ 
plemented  on  the  HP-3000  operating  under  MPE-III. 

7.  External  interfaces 

7.1  System  subroutines 

ATAN 

SQRT 

7.2  Other  programs  called 
None. 


7.  3  External  storage  used 


None. 


8.  Performance  specifications 

8.1  Storage 

Stack:  2  words 
Code:  73  words 

8.2  Execution  time 


8.3  I/O  Load 

There  are  no  I/O  statements. 

8.4  Restrictions 


9.  Method 

The  arcsine  is  calculated  from  a  formula  based  on  [7], 
sec2  y  -  tan2  y  =7.  (1) 

Manipulation  of  this  equation  produces 

tan2 y  =  sec2 y  -  1  =  — - 1  =  — -n  ^ — 

cos2  y  1  -  sin2  y 

If  x  =  sin  y  then  y  =  sin-1  x,  and 


tan2  (sin-1  x)  = 


1  -  x2 


from  which  we  obtain 


sin  ~  1  x  -  tan 


1  Kr^N 


This  gives  the  correct  result  for  x  ^  0.  When 
x<0,  the  relation  [7] 
sin~'  ( -x)  =  -sin-1  x 
is  used. 


It  is  noted  that  when  |x  |  **  1  the  denominator  on  the 
right  side  of  Equation  (2)  is  very  small.  Considerable  loss 
of  precision  due  to  roundoff  error  (in  the  denominator 
itself)  is  expected  in  that  case.  To  improve  the  precision 
of  the  result  as  |x  |  —  1 ,  the  argument  range  is  reduced 
to  [0,  A]  by  the  identity  [9] 

sin-’  x  =  ir/2  -  2  iin~1  ^  ,  (4) 

along  with  Equation  (3)  for  negative  arguments. 

The  sine  of  any  argument  is  restricted  to  the  range 
[-1,  1].  If  x>l  or  x<  - 1,  ASIN  uses 


ASIN  =  k/2 
along  with  Equation  (3). 

ASIN  returns  the  principal  value  of  the  arcsine.  That 
is,  if  y  =  sin  ~ 1  x, 

0  <y  ^  v/2  for  x^O 

-  k/2  ^  y  <0  for  x  <0 

10.  Comments 

None. 


1.  Name 


7.2  Other  programs  called 


MNVAR  None. 


2.  Purpose 

The  purpose  of  MNVAR  is  to  calculate  the  mean  and 
variance  of  a  normalized  univariate  distribution. 

3.  Calling  sequence 

CALL  MNVAR  (HIST,  N,  XMIN,  XMAX, 
XMEAN,  VAR). 

where: 


HIST  =  empirical  density  function  (real) 

N  =  number  of  cells  in  histogram  (integer) 

HIST 

XMIN  =  lower  end  of  variable  range  (real) 

XMAX  =  upper  end  of  variable  range  (real) 

XMEAN  =  mean  of  distribution  (real) 

VAR  =  variance  of  distribution  (real) 


4.  Input— output 

4.1  Input 

There  are  no  input  statements. 

4.2  Output 

There  are  no  output  statements. 

4.3  File  Storage 
None. 

5.  Exits 

There  are  no  nonstandard  exits. 

6.  Usage 

MNVAR  is  written  in  FORTRAN  IV  and  is  presently 
implemented  on  the  HP-3000  operating  under  MPE-III. 

7.  External  interfaces 

7.1  System  subroutines 

None. 


7.3  External  storage  used 
None. 

8.  Performance  specifications 

8.1  Storage 

Stack:  10  words 
Code:  52  words 

8.2  Execution  time 
TBD. 

8.3  I/O  Load 

There  are  no  I/O  statements. 

8.4  Restrictions 
None. 

9.  Method 

The  mean  m  of  a  discrete  probability  distribution  is  the 
first  moment  of  the  distribution, 

m  =  2  x*  P*  (1) 

s 

where  xs  is  one  of  the  values  taken  on  by  the  random 
variable  X  and  ps  is  a  value  of  the  probability  density 
function  of  the  distribution. 

ps  =  Pr{X  =  xs}  (2a) 

£/>,  =  /  (2b) 

s 

The  variance  a 2  is  the  second  moment  about  the  mean 
of  the  distribution, 

-  mf  ps  (3) 
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1.  Name 


7.  External  interfaces 


NORM 

2.  Purpose 

The  purpose  of  NORM  is  to  transform  an  angle  out¬ 
side  the  range  - 180°  <  angle  <  180°.  into  the  com¬ 
parable  angle  within  that  range. 

3.  Calling  sequence 

X  =  NORM  (ANGLE), 
where: 

ANGLE  =  original  angle  (real) 

NORM  =  transformed  angle  (real) 

4.  Input— -output 

4.1  Input 

There  are  no  input  statements. 

4.2  Output 

There  are  no  output  statements. 

4.3  File  Storage 
None. 

5.  Exits 

There  are  no  nonstandard  exits. 

6.  Usage 

NORM  is  written  in  FORTRAN  IV  and  is  presently 
implemented  on  the  HP-3000  operating  under  MPE-III. 


7.1  System  subroutines 
None. 

7.2  Other  programs  called 
None. 

7.3  External  storage  used 
None. 

8.  Performance  specifications 

8.1  Storage 

Stack:  0  words 
Code:  30  words 

8.2  Execution  time 
TBD. 

8.3  I/O  Load 

There  are  no  I/O  statements. 

8.4  Restrictions 
None. 

9.  Method 

If  ANGLE  >  180°  then  360  is  subtracted,  repeatedly 
if  necessary,  until  ANGLE  ^  180®.  Similarly,  if  ANGLE 
^  - 180°  then  360  is  added  as  many  times  as  necessary. 
No  action  is  taken  if  ANGLE  is  already  in  the  correct 
range. 

10.  Comments 
None. 


1.  Name 


7.  External  interfaces 


NRMLZ 

2.  Purpose 

The  purpose  of  NRMLZ  is  to  normalize  a  two- 
dimensional  histogram. 

3.  Calling  sequence 

CALL  NRMLZ  (HIST,  NI,  NJ). 
where: 

HIST  =  histogram  array  (real) 

NI  =  number  of  rows  in  HIST  (integer) 

NJ  =  number  of  columns  in  HIST  (integer) 

4.  Input — output 

4.1  Input 

There  are  no  input  statements. 

4.2  Output 

There  are  no  output  statements. 

4.3  File  Storage 
None. 

5.  Exits 

There  are  no  nonstandard  exits. 

6.  Usage 

NRMLZ  is  written  in  FORTRAN  IV  and  is  presently 
implemented  on  the  HP-3000  operating  under  MPE-III. 


7.1  System  subroutines 
None. 

7.2  Other  programs  called 
None. 

7.3  External  storage  used 
None. 

8.  Performance  specifications 

8.1  Storage 

Stack:  6  words 
Code:  50  words 

8.2  Execution  time 
TBD. 

8.3  I/O  Load 

There  are  no  I/O  statements. 

8.4  Restrictions 
None. 

9.  Method 

All  of  the  elements  of  HIST  are  summed.  The  sum 
is  used  to  normalize  HIST  so  that  the  sum  of  all  the  nor¬ 
malized  elements  =  1. 

10.  Comments 
None. 


1.  Name 


7.2  Other  programs  called 


RMSARY 

2.  Purpose 

The  purpose  of  RMSARY  is  to  find  the  rms  value  of 
a  discrete  univariate  distribution. 

3.  Calling  sequence 

CALL  RMSARY  (ARRAY,  N,  XMIN,  XDELT, 
ARMS). 

where: 

ARRAY  =  unnormalized  univariate  distribution  (real) 

N  =  number  of  cells  in  distribution  (integer) 

XMIN  =  lower  end  of  variable  range  (real) 

XDELT  =  independent  variable  increment  (real) 

ARMS  =  rms  value  (real) 

4.  Input — output 

4.1  Input 

There  are  no  input  statements. 

4.2  Output 

There  are  no  output  statements. 

4.3  File  Storage 
None 

5.  Exits 

There  are  no  nonstandard  exits. 

6.  Usage 

RMSARY  is  written  in  FORTRAN  IV  and  is  present¬ 
ly  implemented  on  the  HP-3000  operating  under  MPE-ID. 

7.  External  interfaces 

7.1  System  subroutines 

None. 


None. 

7.3  External  storage  used 
None. 

8.  Performance  specifications 

8.1  Storage 

Stack:  8  words 
Code:  51  words 

8.2  Execution  time 
TBD. 

8.3  I/O  Load 

There  are  no  I/O  statements. 

8.4  Restrictions 
None. 

9-  Method 

The  rms  value  is  calculated  from 

a  =  [E  x2sf*/N T 

L  s  J 

where  fs  is  a  frequency  distribution  that  satisfies 

£/,  - « 

s 

Here  ARRAY  is  taken  to  represent  {fs}. 

10.  Comments 
None. 
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THE  SUN  SLITTER  ANALYSIS  PROGRAM  IS  NAMED  GLITTER?. 

IT  IS  CONSTRUCTED  FROM  THE  MAIN  PROGRAM  IN  THE  USL  6MATN2RB 
ANO  THE  SUBROUTINES  IN  RL  68UBSRL  AS  FOLLOWS! 


tPREP  GMAIN?RR,GLITTER?»RL*GSUBSRL 
ISAVE  SLITTER2 


THE  SOURCE  FILE  FOR  THE  MAIN  PROGRAM  IS  6MAIN2. 

THE  SOURCE  FILES  FOR  THF  SUBROUTINES  HAVE  THE  SAME  NAMES  AS 
THE  ENTRY  POINTS  IN  GSUH3RL  WITH  THE  FOLLOWING  EXCEPTIONS! 


STEP2A  INSTEAD  OF  STEP? 
STE»5A  INSTEAD  OF  STEPS 
MLACOS  INSTEAO  OF  ACOR 
MLASIN  INSTEAO  OF  ASIN 


THE  SUBROUTINES  ARE  ALSO  IN  THE  USL  FILES  SUBSRB  AND  5U6S2RB , 


GLITTER2  READS  IMAGE  DATA  FROM  LOGICAL  UNIT  1  AND  PRINTS  ON 
UNIT  7.  SO  THE  JOB  SETUP  IS 


!FILE  FTN07i0EV*LPjCCTl 

JFILE  FTN01SF0RMAL  file  DESIGNATOR. OLD 

! RUN  GLITTER? 

ENTER  DATA  In  RESPONSE  TO  PROMPTS  FROM  PROGRAM. 


DATA  ITEMS  ARE  (ALL  ANGLES  TN  DEGREES) ! 


LATITUDE.  LONGITUOF.  OFCLINATION,  TIME 

TIME  IS  GMT.  24-MOUR  SYSTEM,  FORM  MHMM.MM... 

THUS  1810.75  MEANS  16  HR.  10.75  MIN.  PAST 
MIDNIGHT,  OR  6S10J45  PM,  GREENWICH  TIMF. 

AIRCRAFT  HEADING.  ROLL 
WIND  DIRECTION 

CAMERA  FOCAL  LENGTH,  PTCTURE  HEIGHT  (BOTH  SAME  UNITS) 
1ST  ROW,  LAST  RON.  1ST  COLUMN,  LAST  COLUMN 
(IMAGE  AREA  TO  PROCFSS) 

IMAGE  FILE  NAME  (TO  LABEL  PRINTOUT) 


FOLLOWING  AWE  PROGRAM  FILE  LISTINGS  AND  A  SAMPLE  RUN  STREAM. 


nnnnnnnnn 


SCONTROL  MAP. CROSSREF.L ABEL. FTLE*1 

MAIN  PROGRAM  TO  CALCULATE  SEA  SURFACE  SLOPE  HISTOGRAMS! 

1.  alpha-beta  histogram. 

2.  CROSSHINO-UPWINO  histogram. 

MATTHEW  LVBANON,  CSC.  FEBRUARY  15.  1980. 

MODIFIED  FEBRUARY  21.  1880. 

MODIFIED  FEBRUARY  27.  1980. 


PARAMETER  NALPHAsSfc,  NBETA=90 
PARAMETER  NZsai.  7MA*=1.0 

integer  buffer (5i  ?) .  row(2.512),  dalpha,  dbeta 

REAL  LONG,  LAT.  MU.  NU.  H I  ST  ( N  ALPH  A  ,  NEl  E  T  A  )  .  HISTB(NBFTA) 
REAL  HISTCUIN7.N7) ,  SLICE(NZ) 

DOUBLE  PRECISION  COSOMG.  OROLL.  DRROLL.  DOR.  CDR.  SOR. 

♦  CDPHI.  SOPH  T 

CHARACTER*?  FILE 
CHARACTE»«10  NAMFR,  NAME  7 
DATA  NHOW  /3/ 

DATA  DDR  /57. 29577951 508232100/ 

DATA  At,  A?.  81.  B2  /-180.,  180..  0..  90./ 

DATA  NAMES  /"  BETA  " / 

DATA  NAME7  /"  7  "/ 


READ  IN  CONSTANTS. 

WRITE  (8.10001 

1000  FORMAT  ("  ENTER  ALL  ANGLES  IN  DEGREES.*// 

♦  "  ENTER  LATITUDE.  LONGITUDE.  DECLINATION.  TIME  (REAL)  " 
ACCEPT  LAT.  LONG.  DECL.  TIME 

WRITE  (6.1001) 

1001  FORMAT  f"  ENTER  AIRCRAFT  MEAOING.  ROLL  (REALl  "1 
ACCEPT  HOG.  POLL 

WRITE  (6.10051 

1005  FORMAT  ("  ENTER  OIRECTTON  FROM  WHICH  WIND  IS  BLOWING,"/ 

♦  "  MFASUREO  CLOCKWISE  FROM  NORTH  (REAL)  *1 
ACCEPT  WIND 

WRITE  (6,100?) 

1002  FORMAT  ("  ENTER  CAMERA  FOCAL  LENGTH,  PICTURE  HEIGHTswTOTH"/ 

♦  "  (REAL,  BOTH  IN  SAME  UNITS)  ") 


1 


56 

57 


ACCEPT  F.  HEIGHT 
WRITE  (6.10031 

1003  FORMAT  ("  ENTER  FIRST  %  LAST  ROW,  FIRST  &  LAST  COLUMN", 
♦  "  (INTEGER)  ") 

ACCEPT  I R 1 ,  IR2.  TE1.  IF2 
WRITE  (8. 1009) 

100«  FORMAT  ("  ENTER  IMAGE  FILE  NAME  "1 
ACCEPT  FILE 
:  INITIALIZE. 

DU  10  J  =  1 . NBE T A 
HI  STB ( J)  =  0. 

00  10  Isi.NALPHA 
MIST(I.J)  =  0. 

CONTINUE 
00  15  1*1 . N7 
00  15  J*1 .N? 

HISTCU(I.J)  s  0. 
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15  CONTINUE 

i&tci  *  (Iri  ♦  n  /  2 

1SKIP  *  NHOR  ♦  1HECI  -  1 
SCALE  s  HEIGHT  /  SIP. 
n»LPHA  a  360  /  NALPHA 
DhETA  a  90  /  NHETA 
DROLL  «  ROLL 
ORROLL  *  OROLL  /  O0R 
COR  a  COS  (ORROLL 1 
SUR  a  SJN  (ORROLL) 

ZDELT  a  2.  *  7MAX  /  FLOAT  ( N 7 1 

NCNTR  a  (N7  ♦  1)  /  2 

CONVERT  F  TO  SCAN  LINF  UNITS. 

FF  a  F 

F  a  F  /  SCALE 
;  SKIP  TO  START. 

REMIND  I 

OU  20  1*1, I SKIP 

re ao  cn 

20  CONTINUE 

WRITE  (T. 20001  LAT.  LONG.  OECL,  TIME.  HOG.  ROLL.  WIND.  FF .  HEIGHT 

2000  FORMAT  <"l".61X«r0NDTTT0NS'/  .  , 

«,  -OLATITUDE  *»,F8.2.10X"LON6ITUDE  ='.F8.2, 

«  10*"OECLINATION  *" .FT  .2. tOX'TIME  ■".FB.2, 

♦  «  HOURS  (GMT  )  "  / 

+  "OAIRCRAFT  HEADING  * ■ , F8 .2 . 1 0 * 'ROLL  ANGLE  «",F7.2/ 

♦  "Owl NO  DIRECTION  *".F6.2/  .  „ .  . , 

+  "OCAMERA  FOCAL  LENGTH  *" • F4. 1 # 1 0X"P1CTUPE  HEIGHT  «  #Fa.ll 

WRITE  (7.2001)  IRt#  tR2>  I£>»  1£?*  FIl-£  .  .  . 

200!  FORMAT  { "OROWS" , I# , "  -'.I«.'.  COLUMNS'. I#.' 

♦  "OF  IMAGE*. *R) 

WRITE  (7,20021 

2002  FURMAT  ( / / /(.OX "ERROR  MESSAGES'l 
r  CALCULATE  SOLAR  ELEVATION  and  AZIMUTH. 

call  STEP1  (TIME.  LONG#  LAT#  OECL*  PHI#  CDPMI#  SOPH!*  A7#  S300) 

C  CALCULATE  CORRECTION  TO  REFERENCE  ALPHA  TO  WIND. 

DELTA  a  180.  -  (WINO  ♦  *Z) 

C  READ  A  RECORD  (?  SCAN  LINES). 

JO  RE  AO  (1)  BUFFER 

C  UNPACK  SCAN  LINES. 

CALL  GETPIX  (BUFFER.  WOW) 

C 

C  START  LOOP  ON  ROWS. 

00  100  IalRl.IRP 
C  CALCULATE  T. 

y  a  256.5  -  FLOAT  (I) 

C  START  LOOP  ON  COLUMNS. 

00  70  JalEl.IEP 
C  CALCULATE  X. 

X  a  FLOAT  (J)  -  256.5 

C  CALCULATE  MU  ANO  NU. 

CALL  STEP?  (X,  T.  F.  HOG.  AZ.  MU.  NU.  COR.  SDR.  IFLAG) 

C 

IF  (IFLAG  .FO.  >)  GO  TO  J01 
f  CALCULATE  COS  (OMFGA). 

UQ  CALL  STEPS  (TOPHI.  SOPHI,  MU.  NU,  COSOMG) 

C  CALCULATE  BETA. 

CALL  STEP*  (Mil.  SOPHI,  COSOMG.  SETA) 

C  CALCULATE  alpha.  ..... 

CALL  STEPS  (NU.  SETA,  PHI,  COPHI,  COSOMG.  MU,  ALPHA,  S302) 

C 

50  IRON  ■  MOO  (I  *  I.  ?)  ♦  1 

C  INCREMENT  ALPHA-BETA  histogram  cell. 

CALL  STEPS  (ALPHA,  BETA,  ROW(IROw.J).  HIST.  NALPHA. 

*  NBETA ,  OALPHA,  08ETA) 

C  INCREMENT  CROSSWINO-UPWINO  HISTOGRAM  CELL. 

CALL  STEP?  (ALPHA,  BETA,  DELTA,  RON ( IRON, J1 ,  HISTCU. 

♦  N7.  70ELT.  7MAX) 

C  END  LOOP  ON  COLUMNS. 

70  CONTINUE 


29 

JO 

31 

C 

RFAO  NEXT  R£C 090  WHEN  NECESSARY. 

IF  (MOO  fl.  21  .EO.  1)  GO  TO  100 
READ  (l.ENOslOOl  BUFFER 

32 

33 

C 

UNPACK  SCAN  LINES. 

CALL  GETPIX  f BUFFER#  flUw) 

34 

C 

ENO  LOOP  ON  ROWS. 

JS 

100 

CONTINUE 

36 

37 

C 

NORMALIZE  histograms. 

CALL  NWMLZ  (MIST.  NALPHA.  NttETA) 

138 

C 

139 

CALL  NRMLZ  (HISTCU.  N7 .  N 7 1 

140 

C 

FINO  DISTRIBUTION  OVER  BETA. 

141 

00  130  Jsl.NBETA 

142 

DO  130  1=1. NALPHA 

143 

HISTB(J1  *  HISTBf.Il  ♦  HIST  (I.J) 

44 

1  30 

CONTINUE 

145 

C 

FINO  MEAN  AND  VARIANCE  OF  HISTH. 

46 

CALL  MNVAR  (HIST8.  NRFTA.  81.  82.  XMFANB.  VARB) 

147 

C 

FINO  LARGEST  SIGNIFICANT  BETA. 

148 

JMAX  =  0 

149 

JJ  =  N0ETA  ♦  1 

150 

00  150  1=1, NALPHA  * 

151 

00  140  Jsl.NBETA 

152 

K  =  JJ  -  J 

153 

IF  (HISTfl.Kl  .LT.  1.0E-5)  GO  TO  140 

154 

JMAX  Z  MAX  (JMAX.  K) 

155 

GO  TO  150 

156 

140 

CONTINUE 

157 

JMAX  =  MAX  (JMAX.  11 

158 

150 

CONTINUE 

159 

82  s  JMAX  *  OBETA 

160 

C 

PRINT  OUT  CONDITIONS  AND  HISTOGRAMS. 

161 

WRITE  (7.20001  LAT.  LONG.  OECL.  TIME.  HOG.  ROLL.  WIND,  FF ,  HEIGHT 

162 

WRITE  (7.2001;  IR1,  IR2.  IE1.  IE2.  FILE 

163 

WRITE  (7.20041  A 7 

164 

2004 

FORMAT  ("OALPHA  IS  MEASURED  CLOCKWISE  FROM  AN  AXIS 

THAT  ". 

165 

►  "POINTS  TOWARO  THE  SUN."/ 

166 

►  10X-SOIAR  A7IMUTH  *".F7.1,«  COUNTER-CLOCKWISE 

FROM  SOUTH. "1 

167 

WRITE  (7,20031 

/ /65X "BET  A " ) 

168 

2003 

FORMAT  (///50X"N0RMALI7E0  ALPHA-BETA  DISTRIBUTION" 

169 

11  =  NALPHA  /  2  -  ? 

170 

1?  =  11  ♦  6 

171 

no  160  1=1.11 

17? 

write  (  7.20071  (HlSTd.  J)  .  J=l.  JMAXl 

173 

2007 

FORMAT  ( / (6X 1 SF7 .51 1 

174 

160 

CONTINUE 

175 

WRITE  (7,20081  (HISTdltl  .  J)  .  J*1  .JMAX) 

176 

2008 

FORMAT  (3X"A"/(6XI8F7.51 1 

177 

WRITE  (7,2010)  (HlSTdl+2.  J)  .  Jsi.JMAX) 

178 

2010 

FORMAT  (3X"L*/(6X18F7.5) 1 

174 

WRITE  (7.20U)  (HIST  (I1  +  3.J)  »J  =  1  .JMAX) 

180 

201  1 

FORMAT  (3X"P"/(6X18F7.51) 

181 

WRITE  (7,201?)  (HlSTdl"4,J),  Jsi.JMAX) 

18? 

2012 

FORMAT  ( 3X"H«/ (6X1 8F7.5) ) 

183 

WRITE  (7,2008)  (HIST! I1t5. J) . J*1 . JMAX) 

184 

00  170  1=12, NALPHA 

185 

WRITE  (7.20071  (HISTd.  J).J  =  l,  JMAX) 

186 

170 

CONTINUE 

187 

WRITE  (7,2009)  A1.  A2,  OALPHA,  81,  82.  DBETA 

188 

2009 

FORMAT  (//"  ALPHA  RANGES  FHOM»,F6.0."  TO".F6.0. 

1  89 

♦  "  IN  STFPS  OF". 13// 

190 

♦  "  BETA  RANGES  FROm",F6.0« "  TO".F6.0."  IN  STEPS  OF". 13) 

191 

WRITF  (T.2005)  (HISTB(J) . J=1 . JMAX) 

192 

2005 

FORMAT  ("1",55X"0ISTRI8UTI0N  OVER  BETA"// 

193 

♦  (6X18F7.5)) 

194 

XSTEP  «  08ETA 

195 

C 

196 

CALL  HOSPLV  (MI STB ,  NBET A ,  NAMES,  Bl,  XSTEP) 

197 

C 

198 

WRITE  (T, 20061  XMfANB.  V ARB 

199 

2006 

FORMAT  (//10X"MFAN  « • . E 1 5 . 8 . 1 0 X " VAR I ANCE  *".E15.8) 

I* 


200 

WHITE  (7,20131 

20! 

2013  FORMAT  (  "  1  ■ , 46 3 "NORMAL I ZEO  CROSSWIND-UPWIND  DISTRIBUTION"// 

202 

♦  6  43 " UPW I  NO  *  1 

203 

11  *  NON  T9  -  5 

* 

204 

12  «  NCNT9  ♦  5 

205 

DO  100  Isl.Il 

206 

WHITE  (7.20071  (HlSTCU(t.J).Jsl.NZl 

207 

100  CONTINUE 

200 

WHITE  (7.20141 

209 

2014  F09MAT  ,(/33"C»1 

210 

WHITE  (7,20151  f HJSTCtK  1 1 ♦ I »  J) «  J*I » NZ 1 

211 

2015  FORMAT  (■♦■.5310F7.5/(6XlflF7.511 

• 

212 

WRITE  (7,20161 

213 

2016  FORMAT  { / 33 "R" 1 

214 

WRITE  (7,20151  (HISTCII(U+2,Jl.Jcl,N71 

215 

WHITE  (7,20171 

216 

2017  FORMAT  (/33"0"1 

217 

WRITE  (7,20151  (HISTCIKI1W3.J), Jsl.NIl 

210 

WRITE  (7.20101 

219 

2010  FORMAT  (/33*S"1 

• 

220 

WRITE  (7.20151  (HI3TCU(I1*4. J) . Jsl.NZ) 

221 

WRITE  (7.20101 

222 

WHITE  (7.20151  (HISTCU(I1*S.J1,J*1.NZ) 

223 

WRITE  (7,20191 

224 

2019  FORMAT  (/33"W»1 

225 

WRITE  (7.20151  (MISTCU(11»6. J), jsl ,N71 

226 

WRITE  (7.20201 

227 

2020  FORMAT  (/ S3 "  I  *  1 

c 

220 

WRITE  (7.20151  (HTSTCU(Il»7.Jl.Jsl,N71 

229 

WRITE  (7.20211 

230 

2021  FORMAT  ( / 3  3  *N" 1 

231 

WRITE  (7,20151  (HTSTCumwe.Jl.Jsl.N71 

232 

WHITE  (7.20221 

233 

2022  FORMAT  (/33"0"1 

234 

WHITE  (7.20151  (HISTCIJ(T1+9,J1.J»1,NZ1 

235 

00  190  1*12. N7 

• 

236 

WRITE  (7,20071  (HlSTCUdtJ1.Jsl.NZl 

237 

190  CONTINUE 

230 

ZMTN  s  -7MA3 

239 

WHITE  (7,20231  7MIN,  7MA3,  ZDELT 

240 

2023  FORMAT  (//■  BOTH  73  ANO  7Y  RANGE  FROM" , F6 . 2 . »  T0-.F6.2, 

241 

♦  "  IN  STEPS  0F".F6.4) 

242 

C  CROSSWIND  SLICE. 

243 

DO  200  Is(.N7 

• 

24a 

SLICEdl  s  HISTCU(I.NCNTH) 

245 

200  CONTINUE 

246 

WRITE  (7,20241 

247 

2024  FORMAT  ( " 1 » . 553 "CROSSWI NO  0 I S T H IRUT I ON" / 1 

240 

C 

249 

CAUL  HOSPLY  (SLICE,  N7.  NAME Z ,  ZMIN,  ZDELT1 

250 

C 

251 

CALL  RMS ARY  (SLICF,  N7.  ZMIN,  70ELT,  RMS) 

c 

252 

C 

253 

WRITE  (7.20251  RMS 

254 

2025  FORMAT  (//t  03"RMS  SLOPE" . E 1 5. 81 

255 

C  UPWINO  SLICE. 

256 

00  210  J*1.NZ 

257 

SLICE  (J1  s  HISTCU(NCNTR.J) 

250 

210  CONTINUE 

259 

WRITE  (7.20261 

c 

260 

2026  FORMAT  (" 1 ",  56  3 "UPWINO  DISTRIBUTION"/) 

261 

C 

262 

CALL  HPSPLY  (SLICE.  N7.  NAMEZ,  ZMIN,  ZDELT1 

26  3 

C 

264 

CALL  RMSARY  (SLICE,  NZ.  ZMIN,  ZOELT,  RMS1 

265 

C 

266 

WRITE  (7,20251  RMS 

267 

C 

c 

260 

STOP 

269 

C 

270 

C  ERROR  MESSAGES. 

.*■  .■*  .* 


*-*  “J.  V.s 


*,*  \*  *,*  »,*  •,*  *%*  *, 


c 


33 


300  WHITE  {7*3000!  PHI.  *7 

3000  F OHM* T  ("OELFVAT  ION  *".£15.8.10*. 

♦  "A7IMUTH  CAN'T  BE  CALCULATED,  SET  TG".E15.8) 

GO  TO  30 

301  IT  *  ?57  -  Y 

JJ  *  ?57  ♦  * 

WHITE  (7.3010)  II.  JJ.  MU.  NU 

3010  FORMAT  ("OrtAO  MU  CALCULATION  FOR  ROW",  I A .  10*  "COLUMN* ,  1 41 , 

♦  1 0 X "Mu  s".E15.8.!0*"NU  *",£15.8! 

GO  TO  <10. 

30?  IT  *  ?57  -  Y 

JJ  *  ?57  ♦  * 

WRITE  ( 7 ,  30?O )  II.  JJ.  NU.  bETA.  PHI.  COSOMG.  MU.  ALPHA 
30?0  FORMAT  ("ONO  SOLUTION  FOR  ALPHA  EXISTS  FOR  ROW",I«. 

♦  10*"C0LUMN", 14/ 

♦  1 0  X "NU  =".E1S.8.10x"BETA  =  " . E 1 5 . 8 . 1 0 X "PH  I  =».E15.8/ 

♦  10X-COS  (OMEGA)  =".015.8, 10X"MU  =",E15.8/ 

♦  10*" ALPHA  3ET  TO "  , E 1 5 . 8 ) 

GO  TO  50 

C 

END 


v.y.y.v’.v. y.y. 


1  JCONTROL  MAP, CROSSREF. LABEL 

1.1  SilHROUT  I NE  STEP1  (TIME,  LONG.  LAT,  DECL.  PHt.  CDPHI.  SDPHI.  A 

2  C 

3  c  subroutine  to  calculate  local  hour  angle  cuseo  in  subsequent 

a  C  CALCULATIONS! •  SOLAR  ELEVATION,  AND  SOLAR  AZIMUTH. 

5  C 

ft  C  MATTHEW  LYBANON.  CSC.  FEBRUARY  6.  1980. 

7  C 

ft  C  TIME  a  GREENWICH  MEAN  TIME  (MILITARY  FORMAT) 

9  C  LONS  a  LONGITUDE  (♦  FOR  WEST  OF  GREENWICH) 

10  C  LAT  a  LATJTUOE  (♦  FOR  NORTH  OF  EQUATOR) 

11  C  OECL  a  SUN’S  DECLINATION  (♦  FOR  NORTH  OF  EQUATOR) 

12  C  PHI  s  SOLAR  ELEVATTON 

13  C  CDPHI  a  COS  (PHI)  —  00U8LE  PRECISION. 

13.1  C  SDPHI  a  SIN  (PHI)  --  DOUBLE  PRECISION. 

13.2  C  AZ  s  SOLAR  AZIMUTH 

la  C  NOTE!  ALL  ANGLES  IN  ARGUMENT  LIST  ARE  IN  DEGREES. 

15  C  »  *  ERROR  RETURN  TF  AZIMUTH  CAN'T  BE  CALCULATED 

16  C 

17  REAL  LONG.  LAT 

1ft  DOUBLE  precision  dlat.  orlat.  doecl.  oroecl.  oha.  drha, 

19  ♦  NUM.  DENOM.  OOR.  DBPHI,  CDPhI.  SOPHI 

20  DATA  OR.  OOR  /57. 29578.  5 7 . 2 957 7 R5 1 3082 32 1  DO/ 

21  C 

22  C  CALCULATE  LOCAL  HOUR  ANGLE. 

23  I HR  *  TIME  /  100. 

24  XMIN  a  MOD  (TIME.  100.) 

25  C  XTIME  IS  TIME  PAST  MIDNIGHT  IN  MINUTES. 

26  XTIME  s  FLOAT  (60  •  I HR )  ♦  XMIN 

27  C  HA  IS  HOUR  ANGLE. 

28  HA  *  (XTIME  -  720.)  /  4.  -  LONG 

29  C  CALCULATE  SOLAR  ELEVATTON. 

30  C  CONVERT  angles  TO  DOUBLE  PRECISION. 

31  OLAT  a  LAT 

32  DOECL  a  DEC L 

33  OHA  s  HA 

34  C  CONVERT  ANGLES  TO  RADIAN  MEASURE. 

35  ORLAT  s  DLAT  /  DDR 

36  OROECL  *  DOECL  /  OOR 

37  ORHA  a  OHA  /  OOR 

38  C 

39  Z  *  SIN  (ORLAT)  *  SIN  (OROECL)  ♦  COS  (ORLAT)  * 

40  ♦  COS  (ORDECL)  •  COS  (ORHA) 

41  C  PHI  IS  SOLAR  ELEVATION  (RAO  I  AN  MEASURE). 

42  PHI  5  ASIN  (Z) 

43  C  CALCULATE  SOLAR  AZIMUTH. 

44  C  CONVERT  PHI  TO  DOUBLE  PRECISION. 

45  ORPHI  a  PHI 

46  CDPHI  a  COS  (ORPHI) 

46.1  SDPHI  s  SIN  (ORPHI) 

46.2  NUM  a  -COS  (OROECL)  *  SIN  (ORHA) 

47  DENOM  a  CDPHI 

4ft  IF  (ABS  (NUM)  .LE.  ABS  (OENOM))  GO  TO  10 

49  IF  (OENOM  .NE.  0.000)  GO  TO  10 

50  C  AZIMUTH  CANNOT  HE  CALCULATED  (SUN  IS  PROBABLY  AT  ZENITH). 

51  C  SET  TO  ZERO  AND  TAKE  FRROR  EXIT. 

52  A7  a  0. 

53  PHI  a  PHI  *  OR 

54  RETURN  1 

55  C  AZIMUTH  CAN  BE  CALCULATED. 

56  10  Z  a  NUM  /  OENOM 

57  C  hi  IS  SOLAR  AZIMUTH. 

58  AZ  a  ASIN  (7) 

59  c  CONVERT  PHI  AND  AZ  TO  DEGREES. 

60  PHI  a  PHI  •  OR 

61  AZ  a  AZ  *  OR 

62  C 

65  RETURN 

64  C 
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21 
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20 

29 
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31 

32 

33 
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35 

36 

37 
30 

39 

40 

41 
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SCONTHOL  MAP.CROSSREF.LABEl, 

SUBROUTINE  STEP?  (X,  Y.  F.  HOG.  *7.  Mu.  NIJ.  CUR.  SDR.  IFLAG) 
C 

C  SUBROUTINE  TO  CALCULATE  IMAGE  ANGULAR  COORDINATES. 

C  VERSION  INCLUOING  AIRCRAFT  ROLL  COMPENSATION. 

C 

C  MATTHEW  LYBANON.  CSC.  FEBRUARY  6,  1900. 

C  MOOIFIEO  FEBRUARY  15.  1980. 

C 

C  X 

C  Y 

C  F 

C  NOTE! 

C  HOG 

C  A  7 

C  MU 

C  NO 

C 

C  COR 

C  SOR 

C  IFLAG 

C 
C 
C 

C  note: 

C 

REAL  MU.  Nil.  NORM 

DOUBLE  PRECISION  ORMU.  07.  OX,  OY.  OF,  DDR. 

♦  0.  NUM.  OENOM,  COR.  SOR 

OATA  DOR  /57. 29577951308232100/ 

C 

IFLAG  a  0 

C  CONVERT  X.  Y,  ANO  F  TO  DOUBLE  PRECISION. 

OX  *  X 
OY  s  Y 
OF  s  F 
C 

IF  (SDR  .NE.  0.000)  GO  TO  30 
C  7ER0  ROLL  CALCULATION. 

C  CALCULATE  MU. 

07  a  SORT  (OX  •  OX  ♦  OY  •  OY)  /  OF 
ORMU  a  ATAN  (07) 

C  CONVERT  MU  TO  DEGREES. 

M|j  a  ORMU  •  DOR 
C  CALCULATE  NU. 

IF  (X  .NE.  0.1  GO  TO  10 
C  SPECIAL  CASE.  CAN'T  DIVIDE  BY  0. 

IF  (Y  .GE.  0.1  7  a  90. 

IF  (Y  .LT.  0.)  7  =  -90. 

GO  TU  20 

C  NORMAL  CALCULATION. 

10  02  a  ATAN  (OY,  OX) 

C  CONVERT  7  TU  OEGRFES . 

2  a  07  »  DDR 
C 

20 
C 


=  COORDINATE  ALONG  AXIS  POINTING  TOWARD  STARBOARD  WING, 
a  COORDINATE  ALONG  AXIS  POINTING  TOWARD  AIRCRAFT  NOSE, 
a  CAMERA  LENS  FOCAL  LENGTH  (ASSUMED  POSITIVE). 

X,  Y.  AND  F  HAVE  THE  SAME  UNITS, 
a  AIRCRAFT  HEADING, 
a  SOLAR  A7IMUTH. 

a  ANGLE  BETWEEN  REFLECTED  RAY  AND  VERTICAL, 
a  ANGLE  OF  IMAGE  POINT  MEASURED  CLOCKWISE  FROM  AN  AXIS 
POINTING  AWAY  FROM  SUN. 

a  COSINE  OF  AIRCRAFT  ROLL  ANGLE  (DOUBLE  PRECISION), 
a  SINE  OF  AIRCRAFT  ROLL  ANGLE  (DOUBLE  PRECISION), 
a  FLAG  THAT  INDICATES  ERROR  IN  MU  CALCULATION 
FOR  N0N7ERO  ROLL, 
a  0  --  NO  ERROR . 
a  1  —  ERROR. 

ALL  ANGLES  IN  ARGUMENT  LIST  ARE  IN  DEGREES. 


N(l  a  NORM  (A7  ♦  HDQ  ♦  RO 


7) 


V.' 


■■t'.MTyi  ^'■A'A'ivrrA  tt  ffr  i'  r-  u  ij 


1 1 


-rr~i*—r:v 


59 

RETURN 

60 

C 

nonzero  roll  calculation. 

61 

C 

CALCULATE  mu. 

62 

30 

0  3  COR  -  OX  *  SDR  /  OF 

63 

NUM  3  OX  *  OX  ♦  DT  *  07  ♦ 

64 

♦  Of  *  SOW  *  (OF  *  SOW 

65 

IF  (NUM  .LT.  0.000)  IFLAG  3  1 

66 

NUM  3  SORT  (ASS  (NUM) 1 

67 

0EN0M  3  OF  *  0 

68 

0RMU  3  ATAN  (NUM.  OENOM) 

69 

C 

CONVERT  MU  TO  DEGREES. 

70 

MU  *  DRMU  •  00R 

71 

C 

calculate  nu. 

72 

DENOM  s  OF  *  SOW  +  DX  *  COR 

73 

IF  (OENOM  ,NE.  0.000)  GO  TO  40 

74 

C 

SPECIAL  CASE.  CAN'T  OIVIDE  BY 

75 

IF  (7  .GE.  0.)  2  3  90. 

76 

IF  (7  .LT.  0.)  2  *  -90. 

77 

GO  TO  50 

78 

C 

NORMAL  CALCULATION.. 

79 

40 

0?  3  ATAN  (07.  DENOM) 

80 

C 

CONVERT  2  TO  DEGREES. 

81 

2  3  02  *  DOR 

82 

C 

83 

50 

NU  3  NORM  (A?  ♦  HOG  ♦  90.  •  2) 

84 

C 

85 

RETURN 

86 

C 

87 

END 

o*  *  (cow  ♦  on 


C 


t 


C 
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SCONTROL  MAP.CROSSREF.L ABEL 

SUBROUTINE  STEP3  fCOPHI.  SOPHI,  MU,  NU.  CQSOMG) 

SUBROUTINE  TO  CALCULATE  COS  (UMEGA),  WHERE  OMEGA  IS  THE 
ANGLE  OF  INCIDENCE  or  reflection. 

Matthew  lybanon,  esc.  February  b»  ibbo. 

CDPHI  a  COSINE  OF  SOLAR  ELEVATION  -•  DOUBLE  PRECISION. 
SOPHI  a  SINE  OF  SOLAR  ELEVATION  --  DOUBLE  PRECISION. 

Mil  a  ANGLE  BETWEEN  REFLECTEO  RAY  ANO  VERTICAL. 

NU  =  ANGLE  OF  IMAGE  POINT  MEASURED  CLOCKWISE  FROM  AN 

AXIS  THAT  POINTS  AWAY  FROM  SUN. 

NOTES  ALL  ANGLES  IN  ARGUMENT  LIST  ARE  IN  DEGREES. 

COSOMG  a  COS  ( OMEGA )  --  OOUBLE  PRECISION. 

REAL  MU,  NU 

OOUBLE  PRECISION  COPHT.  SOPHI,  DMU,  DRMO.  DNU.  DRNU.  07. 
♦  DDR,  COSOMG 

DATA  DOR  /57.7B5T7B513082321D0/ 

C  CONVERT  ANGLES  TO  DOUBLE  PRECISION. 

OMU  a  MU 
ONU  a  NU 

C  CONVERT  ANGLES  TO  RADIAN  MEASURE. 

DRMU  =  OMU  /  DDR 
ORNU  a  ONU  /  DOR 
C 

07  a  -COPHI  *  SIN  (0RMU1  *  COS  (DRNU)  ♦ 

.  SOPHI  *  COS  (DRMU) 

C  MAKE  SURE. 

07  a  MAX  (DZ ,  -1 .000) 

DZ  a  MIN  (D? •  t  .000) 

COSOMG  a  SORT  (0.500  *  (OZ  ♦  I. 000)) 


C 


RETURN 


to 


i 

i 

3 

tt 

5 
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6 
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10 
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12 

13 

la 

15 

16 
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f CONTROL  MAP, CROSSREF. LABEL 

SUBROUTINE  STEPa  fMU.  SOPHl#  COSOMG.  BETA) 
SUBROUTINE  TO  CALCULATE  SEA  SURFACE  TTLT. 


10 


MA  T  THEW  L YBANON,  CSC.  FEBRUARY  6.  1980. 

*  ANGLE  BETWEEN  REFLECTED  RAT  AND  VERTICAL. 


MU 

SOPHl 

COSOMG 

BETA 
NOTE  J 


*  SINE  OF  SOLAR  ELEVATION  —  DOUBLE  PRECISION. 

*  COSINE  OF  ANGLE  OF  INCIOENCE  OR  REFLECTION 

—  DOUBLE  PRECISION, 
s  SEA  SURFACE  TILT. 

ALL  ANGLES  IN  ARGUMENT  LIST  ARE  IN  DEGREES. 


REAL  MU 

DOUBLE  PRECISION  OMU.  DRMU.  SDPHI,  COSOMG.  DOR 
DATA  OR.  OOR  /57. 29578,  57.29577951308232100/ 

IF  (COSOMG  . NE .  0.0001  GO  TO  10 

GLANCING  INCIOENCE  ••  NO  SULUTION  FOR  BETA. 

ASSIGN  ARBITRARY  VALUE  ANO  RETURN. 

BETA  s  0.0 


RETURN 
CONVERT  MU 
OMU  *  MU 
CONVERT  MU 
DRMU  *  OMU 


TO  00U8LE  PRECISION, 

TO  RADIAN  MEASURE. 

/  OOR 


Z  a  (COS  (0RMU1  ♦  SOPHl)  /  (2.000  *  C0S0MG1 
BETA  9  ACOS  (Z) 

CONVERT  BETA  TO  DEGREES. 

BETA  a  BETA  •  OR 

RETURN 


END 


« 


« 
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IT 

18 
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28 
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30 
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SCONTROL  MAP, CROSSREF, LABEL 

SUBROUTINE  STEPS  (NU.  BFTA,  PHI.  CDPHI,  COSOMG.  MU.  ALPHA.  *) 
C 

C  SUBROUTINE  TO  CALCULATE  LOCAL  AZIMUTH  OF  ASCENT 

C  OF  SEA  SURFACE. 

C  MODIFICATION  1.  SOLUTION  OF  KINSMAN'S  2ND  GRID  EQUATION 

C  USING  CRITERION  TO  RESOLVE  AMBIGUITY. 

C 

C  MATTHEW  LY8ANON,  CSC.  FEBRUARY  7,  1980. 

C  MOUIFIEO  FEBRUARY  13.  1980. 

C 

C  NU 

c 

C  BETA 

C  PHI 

C  COPHI 

C  COSOMG 

C 

C  Mli 

C  ALPHA 

C  NOTES 

C  • 

C 

REAL  mu.  nu 

DOUBLE  PRECISION  ONU.  DRNU .  ObETA.  ORBETA.  CDPHI. 

♦  COSOMG.  ORALP.  OTEST.  A.  8,  02.  NUM.  DENOM,  DDR. 

♦  DMU.  DRMU 

OATA  DR.  DOR  /57. 29575,  57 . 295 7 7951 308232 1 00/ 

C 

C  CONVERT  ANGLES  TO  OOUBLE  PRECISION. 

DNU  a  NU 
OBETA  x  BETA 
DMU  a  MU 

C  CONVERT  ANGLES  TO  RADIAN  MEASURE. 

ORNU  *  ONU  /  DOR 
ORBETA  *  OBETA  /  OOR 
ORMU  x  DMU  /  DOR 
C  CHECK  FOR  SPECIAL  CASES. 

IF  (COPHI  .NE.  0.0001  GO  TO  10 
C  SUN  IS  AT  ZENITH.  SO  ORIENTATION  OF  COORDINATES  IS 

C  UNDEFINED. 

ALPHA  s  NU 

C 

RETURN 
C 

10  IF  (SIN  (DRBETA1  .NF.  0.000)  GO  TO  20 
C  NO  WAVE  TILT.  SO  ALPHA  IS  NOT  OEFINEO. 

ALPHA  s  0. 

C 

RETURN 
C 

20  IF  (COSOMG  .NE.  0.0001  GO  TO  30 
C  GLANCING  INCIDENCE.  NO  SOLUTION  FOR  ALPHA. 

ALPHA  x  0. 

C 


x  ANGLE  OF  IMAGF  POINT  MEASURED  CLOCKWISE  FROM  AN 
AXIS  THAT  POINTS  AWAY  FROM  THE  SUN. 
s  SEA  SURFACE  TILT, 
x  SOLAR  ELEVATION, 
x  COS  (PHI1  -•  OOUBLE  PRECISION. 

=  cosine  of  angle  uf  incidence  or  reflection 
--  double  precision. 

X  ANGLE  BETWEEN  REFLECTED  RAY  ANO  VERTICAL, 
x  SOLUTION  OF  SECOND  GRIO  RELATION. 

ALL  ANGLES  IN  ARGUMENT  LIST  ARE  IN  DEGREES, 
x  ERROR  RETURN  IN  CASE  ALPHA  CAN'T  BE  CALCULATED. 


C 


RETURN 


IF  (SIN  (0RMU)  .NE.  0.000)  60  TO  «0 
REFLECTED  BEAM  IS  IN  ZENITH  DIRECTION. 

ALPHA  3  ISO. 

RETURN 

IF  (SIN  (ORNU)  . NE.  0.000)  60  TO  70 
IMAGE  POINT  IS  ON  Y-AXIS. 

IF  (NU.  .NE.  ISO.)  GO  TO  50 
ALPHA  3  ISO. 

RETURN 
NU  3  0. 

Z  s  COSOMG 

OMEGA  s  ACOS  (Z)  *  OR 

IF  (ABS  (PHI  4  OMEGA  •  BETA)  .LE.  AbS  (PHI  4  OMEGA  4  BETA)) 
►  GO  TO  60 
ALPHA  s  ISO. 

RETURN 


ALPHA  s  0. 

RETURN 

CALCULATE  ALPHA. 

A  x  1.00  /  TAN  (ORNU) 

B  3  COPHI  /  (2.0  «  SIN  (0R8ETA) 
TEST  FOR  EXISTENCE  OF  SOLUTION. 
OZ  3  A  *  A  4  1.000  -  B  *  B 

IF  (02  .GE.  0.000)  60  TO  SO 
NO  SOLUTION. 

ALPHA  s  ATAN  (A  •  B,  -B)  4  DOR 


RETURN  1 

SOLUTION  EXISTS. 

07  s  SORT  (OZ) 

OTEST  3  SIN  (ORNU)  4  (SIN  (DRM U)  - 
NUM  3  A  •  6  4  SIGN  (07.  OTEST) 

OENOM  3  .6  4  A  4  SIGN  (OZ.  OTEST) 

ORALP  3  ATAN  (NUM.  OENOM) 

CONVERT  ALPHA  TO  DEGREES. 

ALPHA  3  ORALP  4  OOR 


COSOMG) 


COS  (ORNU)  4  COPHI) 


RETURN 


v> 
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SCONTROL  M4P.CR0SSREF .LABEL 

•SUBROUTINE  STEPS  f  ALPHA ,  SETA.  POINT,  HIST,  NALPHA ,  NBETA 
♦  OALPHA,  OBETA) 

SUBROUTINE  TO  INCREMENT  ALPHA-BETA  HISTOGRAM. 

MATTHEW  LY8AN0N.  CSC.  FEBRUARY  15,  1BB0. 

ALPHA-  *  AZIMUTH  OF  ASCENT  OF  SEA  SURFACE. 

BETA  x  TILT  OF  SEA  SURFACE. 

POINT  x  IMAGE  POINT  INTENSITY  VALUE. 

HIST  x  HISTOGRAM  ARRAY. 

NALPHA  X  ROW  DIMENSION  OF  HIST. 

NBETA  x  COLUMN  DIMENSION  OF  HIST. 

OALPHA  x  WIDTH  OF  AN  ALPHA  BIN. 

DBETA  x  WIDTH  OF  A  BETA  BIN. 

INTEGER  POINT.  OALPHA.  OBETA 
REAL  HISTINALPHA.NBETA1 
C 

I  X  NALPHA  -  IF  I X  I1SO.  -  ALPHA)  /  OALPHA 
J  x  NBETA  -  I F I X  f RO.  -  BETA)  /  OBETA 
J  x  MAX  U.  1) 

HIST(I.J)  x  HIST(I.J)  ♦  POINT 
C 

RETURN 

C 

END 


[  ------ 

* 

i 

SCONTROL  MAP, CROSSREF. LABEL  I 

2 

SUBROUTINE  STEP7  (ALPHA.  BETA.  OELTA,  POINT.  HIST.  j 

3 

♦  N.  ZOELT .  ZMAX)  i. 

• 

u 

5 

c 

c 

SUBROUTINE  TO  INCREMENT  CROSSwIND-UPwIND  HISTOGRAM. 

6 

7 

c 

c 

MATTHEW  LYB ANON.  CSC.  FEBRUARY  26.  1980. 

a 

4 

c 

c 

ALPHA  *  AZIMUTH  OF  ASCENT  OF  SEA  SURFACE. 

1  0 

c 

BETA  a  TILT  OF  SEA  SURFACE. 

1 1 

c 

DELTA  a  CORRECTION  TO  REFERENCE  ALPHA  TO  (UPWIND) 

• 

12 

c 

WIND  OIRECTION. 

1  3 

c 

POINT  a  IMAGE  POINT  INTENSITY  VALUE. 

1  4 

c 

HIST  a  HISTOGRAM  ARRAY. 

IS 

c 

N  a  ROW  AND  COLUMN  DIMENSION  OF  HIST  (SHOULD  BE  ODD 

16 

c 

SO  ORIGIN  IS  AT  CENTER  OF  A  BIN). 

i  r 

c 

ZOELT  a  MIOTH  OF  A  BIN  (EITHER  AXIS). 

is 

c 

ZMAX  a  UPPER  LIMIT  FOR  EACH  AXIS. 

• 

19 

c 

20 

INTEGER  POINT 

21 

REAL  HIST(N.N) 

22 

DATA  OR  /57. 29578/ 

23 

c 

24 

c 

CORRECT  ALPHA  TO  WIND  DIRECTION. 

23 

AR  a  ALPHA  ♦  OELTA 

c 

26 

c 

CONVERT  ALPHA  ANO  BETA  TO  RADIANS. 

27 

AR  a  AR  /  OR 

28 

8R  a  BETA  /  OR 

29 

c 

calculate  wave  SLOPE. 

30 

SLOPE  a  TAN  (BR1 

31 

c 

CALCULATE  projections  on  histogram  AXES. 

32 

ZX  a  SLOPE  •  SIN  (AR) 

• 

33 

ZY  a  SLOPE  *  COS  (AR) 

34 

c 

CALCULATE  INOICES  OF  HISTOGRAM  CELL. 

35 

I  a  N  -  IFIX  ((ZMAX  -  ZX)  /  ZOELT) 

36 

I  a  MAX  (I.  1) 

37 

I  a  MIN  (I.  N) 

36 

J  a  N  -  IFIX  ((ZMAX  -  ZY)  /  ZOELT) 

39 

J  3  MAX  (J.  1) 

• 

40 

J  a  MIN  (J.  N) 

41 

c 

INCREMENT  HISTOGRAM  CELL. 

42 

HIST ( I . J)  a  HIST ( I  *  J)  ♦  FLOAT  (POINT) 

43 

c 

44 

RETURN 

45 

c 

46 

ENO 

c 

c 

c 

43 

c 

’**•  /'*.  •* 

/-/c  •  •  -.vo  •/ .  *  •  •  *  •  •  •  •  a*  o.o  *  *  v*  - "  *  *  •  •  *  ’  o  •  •  ‘V*  *•'*/» /a’  v\*.vO  o 

2 

SIIBROUT I  NE  GETPI*  f  8l 

3 

C 

4 

c 

SUBROUTINE  TO  UNPACK 

5 

c 

b 

c 

MATTHEW  LY8AN0N.  CSC 

7 

c 

8 

c 

BUFFER  *  READ  BUFFER 

4 

c 

ROW  *  2-ROW  ARRAY 

10 

c 

11 

INTEGER  BUFFER f 51 ?) # 

12 

c 

13 

IOFSET  =  1 

14 

00  100  IROWsl ,? 

15 

00  50  ICOLsl.512.2 

lf> 

ROW ( I ROW ,  ICOL 1  s 

17 

ROW ( I  ROW* ICOL+1 ) 

18 

50 

CONTINUE 

19 

IOFSET  *  513 

20 

100 

CONTINUE 

21 

c 

22 

RETURN 

23 

c 

24 

ENO 

© 


e 


! 


44 


C 


SCONTROL  MAP.CROSSREF.t  ArtEI. 

SliBRQUT  INE  HOSPLY  (HIST.  N,  NAME.  XO.  XSTEP) 

C 

c  SlJBROUT  I NE  TO  DISPLAY  A  UNIVARIATE  HI  STOGRAM. 

C 

C  MATTHEW  L YBANON .  CSC.  FEBRUARY  22.  1980. 

C 

C  HIST  *  HISTOGRAM  ARRAY. 

C  N  *  NUMBER  OF  ENTRIES  IN  HIST. 

C  NAME  *  NAME  OF  INDEPENDENT  VARIABLE. 

C  XO  *  LOWER  ENO  OF  VARIABLE  RANGE. 

C  XSTEP  *  WIDTH  OF  A  BIN. 

C 

REAL  HISTCNI 
CHARACTER*10  NAME 
CHARACTER*!  STAR 
DATA  STAR  /•*■/ 
c 

C  FIND  MAXIMUM-OCCUPANC Y  CELL  AND  CALCULATE  SCALE  FACTOR. 

HMAX  *  HIST(l) 

DO  10  1*2. N 

HMAX  *  MAX  (HMAX.  HIST(I)) 

10  CONTINUE 

FCTR  *  100.  /  HMAX 
C  DISPLAY  HISTOGRAM. 

WRITE  (7,1 000)  NAME 

1000  FORMAT  ("0".A10."  OCCUPANCY"/) 

X  s  XO 

WRITE  (7.1001)  X 

1001  FORMAT  ( 1 XF9 . 3 ) 

DO  40  1*1. N 

NSTAR  3  FCTR  *  HIST(I)  ♦  O.S 

IF  (NSTAR  .LE.  0)  GO  TO  20 

WRITE  (7.100?)  HIST ( I ) .  (STAR, K*I .NSTAR) 

1002  FORMAT  (11XF9.7.IX100AI) 

GO  TO  30 

20  WRITE  (7.100?)  HIST ( I ) 

30  X  s  X  ♦  XSTEP 

WRITE  (7.1001)  X 
40  CONTINUE 
C 

RETURN 


SCONTROL  mar.CWOSSREF. label 
FUNCTION  ACOS  (X) 

c 

C  CALCULATES  THE  INVERSE  COSINE  OF  X 

C 

C  MATTHEW  LYBANON,  CSC.  FEBRUARY  H. 

C 

AC03  a  1.5707R6  -  ASIN  (XI 
C 


RETURN 


c 


scontrol  map. crossref. label 


2 

FUNCTION  A3IN  1X1 

5 

C 

4 

C 

CALCULATES  THE  inverse  SINE  OF  X. 

5 

c 

b 

c 

MATTHEW  LVBANON.  CSC.  FEBRUARY  11. 

7 

c 

8 

IF  (ABS  (X)  .GE.  1.01  0(1  TO  20 

9 

c 

REDUCE  ARGUMENT  RANGE  TO  (0,  0.5). 

10 

IF  (ABS  (XI  .LE.  0.51  GO  TO  10 

11 

c 

ABS  (X)  .GT.  0.5  BUT  .LT.  1.0. 

12 

Y  a  (1.  -  ABS  (XII  /  2. 

1  3 

ASIN  a  1.57079b  -  ?.  *  ATAN  (SORT 

14 

GO  TO  30 

15 

c 

ABS  (X)  .LE.  0.5. 

lb 

10 

ASIN  a  ATAN  (SORT  (X  *  X  /  (1.  - 

17 

GO  TO  30 

18 

c 

ASS  (X)  a  t.O. 

19 

20 

ASIN  a  1.57079b 

20 

C 

TAKE  CARE  OF  NEGATIVE  ARGUMENTS. 

21 

30 

IF  (X  .LT.  0.1  ASIN  a  -ASIN 

22 

C 

23 

RETURN 

24 

C 

25 


END 


» 


» 


I 


> 


I 


1980. 
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1  SCONTRUL  MAP.CROSSREF.LABFL 

2  SUBROUTINE  MN  VAR  fHIST.  N.  XMIN,  XMAx,  XMEAN.  VAR) 

1  C 

4  C  SUBROUTINE  TO  CALCULATE  THE  MEAN  AND  VARIANCE  OF  A 

5  C  NORMALIZED  UNIVARIATE  DISTRIBUTION. 

6  C 

7  C  MATTHEW  LYRANON.  CSC,  FEBRUARY  21,  1980. 

8  C 

9  C  HIST  *  EMPIRICAL  DENSITY  FUNCTION. 


10 

C 

N  s  NUMBER  OF  ELEMFNTS  in  HIST. 

11 

C 

XMIN  S  LOWER  END  OF  VARIABLE  RANGE. 

12 

c 

XMAX  s  UPPER  END  OF  VARIABLE  RANKE. 

1 1 

c 

XMEAN  s  MEAN  OF  DISTRIBUTION. 

14 

c 

VAR  s  VARIANCE  OF  DISTRIBUTION. 

15 

c 

18 

REAL  HIST(N) 

17 

c 

18 

8IIM1  =  0. 

19 

SUM2  r  0. 

20 

OX  s  fXMAX  -  XMIN)  /  FLOAT  IN) 

21 

X  =  XMIN  ♦  0.5  *  OX 

22 

DO  10  I s i ( N 

23 

SUM 1  3  SUM1  ♦  HIST(I)  *  X 

24 

SUM2  s  SUM2  ♦  HISTII)  A  X**2 

25 

X  s  X  ♦  OX 

2b 

10 

CONTINUE 

27 

XMEAN  3  SUM  1 

28 

VAR  s  SUM2  -  SUM1 **2 

29 

c 

10 

RETURN 

11 

c 

12 

ENO 

t 

? 

5 


a 

5 

6 

7 

8 
9 

10 
1 1 
12 
1? 
la 
15 
18 

17 

18 


SCONTROL  MAP. CROSSREF. LABEL 

REAL  FUNCTION  NORM  (ANRI  El 
C 

c  TRANSFORMS  angle  to  THE  RANGt  f-180,  1801. 
C 

C  MATTHEW  L YB ANON .  CSC.  FF6RUARY  7.  1980. 

C 

NORM  X  angle 

10  IF  (NORM  .LE.  180.01  GO  TO  20 
NORM  s  NORM  -  580.0 
GO  TO  10 
C 

20  IF  (NORM  .GT.  -180.01  RETURN 
C 

NORM  =  NORM  ♦  580.0 
GO  TO  20 
C 

ENO 


« 


* 


# 


49 


^n-Ln-.s.- 


nnnnnnnnn 


SCONTROL  MAP.CR0S8RFF.LABEI 

SUBROUTINE  NRML7  (HIST  •  M.  NJl 

SUBROUTINE  TO  NORMAL  T 7F  A  2-DIMENSIONAL  HISTOGRAM 

MATTHEW  LYBANON.  CSC.  FEBRUARY  2*.  1«»S0. 

HT ST  s  HISTOGRAM  ARRAY. 

NT  S  ROW  DIMENSION  OF  HIST. 

N.f  *  COLUMN  DIMENSION  OF  HIST. 

REAL  histcni.njy 

SUM  s  0. 

00  10  Isl.NI 
DO  10  Jsl.NJ 

SUM  s  SUM  ♦  HISTfl.J) 

10  CONTINUE 

00  20  Isl.NI 
DO  ?0  Jsl.NJ 

HIST ( I . J)  s  HIST(I.J)  /  SUM 
20  CONTINUE 
C 


RETURN 


SCONTMm  rtAP.CROSSREF.L ABEl 

SUBROUTINE  RMSARY  f  ARRAY  >  N,  XMIN.  XDFLT.  ARMS! 

C 

C  SUBROUTINE  TO  FIND  THF  RMS  VALUE  OF  AN  ARRAY. 

C 

C  MATTHEM  LVBANON.  CSC.  FEBRUARY  SB.  1BB0. 

C 

C  ARRAY  a  UNIVARIATE  DISTRIBUTION  f  MAY  NOT  BE  NORMAL I7FDI 

C  N  ■  NUMBER  OF  ELEMENTS  IN  ARRAY. 

C  XMIN  a  LOmER  END  OF  VARIABLE  RANGE. 

C  XOELT  a  STEP  SITE. 

C  ARMS  a  RMS  VALUE  OF  ARRAY. 

C 

REAL  ARRAY  INI 
C 

D  a  ARRAY (IT 

X  a  XMTN  ♦  O.S  a  XOFLT 
SUM  a  ARRAY  f 1 )  *  X**E 
DO  10  Ia?.N 

0  a  0  ♦  ARRAY  III 

X  a  X  ♦  XOELT 

SUM  a  SUM  ♦  ARRAY  IT)  *  X**2 

10  CONTTNUE 
C 

ARMS  a  SORT  T SUM  /  0) 

C 

RETURN 


•  *  •  .  • «T  . * V ‘  •  ».y j*  ,t.  •  .  ' 

'  •  '  i  k  /  *  •  */  *•*»’»  *  /  ".*  *, 


_  ■  —  •  a  .  .  *  *  a  a  a  .  •  -  • 


SAMPLE  RUN  STREAM 


! JOB  GLINT. C335 
!  PI  LE  PTN07; DEV-LP; CCTL 
! PI LE  PTNOl- 10059711 . J01 33524 . SI 01 .OLD 
♦RUN  GL1TTBR2 
30 . 84 . 71 . 92 , 6 . 66 , 1606 . 39 
0.  .0. 

100. 

3.  ,4.5 
I  .512,1,512 
10059711 

! P I LE  PTN01 = I 0059582 . J01 33524 . S101 .OLD 
•RUN  GLITTER2 
30.34,72.00,6.66,1704.44 

90. ,-23.66 
100. 

3 .  ,4.5 

1.512.1.512 
1 0059582 

! PILE  PTN01-I0059614.J0133524.S101, OLD 
♦RUN  GLITTER2 
30.34.72.00,6.66,1704.46 
90.  ,-23.66 
100. 

9  4  5 

1 .512.1.512 
10059614 

! PI LE  FTN01-I 0059631 . J01 33524 . S101 ,OLD 
•RUN  GLITTER2 
30.34,72.00,6.66,1704.49 

90. , -23.66 
100. 

3. . 4.5 

1 .512.1.512 
10059631 

!  EOJ 
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Appendix  A:  Solution  for  local  azimuth  of 
sea-surface  ascent 


A1.  Introduction 

This  appendix  discusses  the  solution  for  a,  the  local 
azimuth  of  sea-surface  ascent.  The  equation  for  a  is  suffi¬ 
ciently  complex  that  the  properties  of  the  solution  are  not 
obvious.  One  property  that  will  be  proved  is  that  (when 
there  is  a  solution  at  all)  there  are  two  roots  in  general. 
This  poses  a  problem  in  solving  for  the  (a,  j8)  that  cor¬ 
responds  to  each  pair  of  values  in  the  film  plane  (ji ,  v). 
It  is  to  be  expected  that  further  consideration  of  the  physics 
of  the  situation  will  be  required  to  resolve  the  dilemma. 

A2.  Analysis 

The  azimuth  angle  a  is  a  solution  of  the  equation 

cotv  =  cola  +  lA  csca  csc/3  seco:  cos<f>,  (Al) 

where  v  gives  the  angle  that  the  image  point  makes  with 
the  film  plane  y-axis,  /3  is  the  inclination  of  the  sea  sur¬ 
face,  0)  is  the  angle  of  incidence  or  reflection,  and  <f>  is 
the  solar  elevation.  The  sign  of  the  last  term  is  the  op¬ 
posite  of  that  given  in  Reference  2  (second  of  Equations 
(2)),  but  a  careful  check  of  the  derivation  gives  the  result 
shown  above. 

Equation  (Al)  can  be  simplified  by  multiplying  through 


by  sina.  This  gives 

sina  cotv  =  cosa  +  Vi  csc@  seco:  cos<t>,  (A2) 
which  can  be  written 

A  sina  -  cosa  +  B,  (A3) 

where 

A  =  cotv,  (A4a) 

B  =  l/t  cscfi  seco:  cos<j>.  (A4b) 


The  properties  of  the  solution  of  Equation  (A3)  can  be 
derived  by  writing  the  equation  in  a  different  form.  Make 
the  change  of  variables 

x  =  cosa  (A5a) 

y  =  sina  (A5b) 


With  these  substitutions  Equation  (A3)  becomes 


y  =  x/A  +  B/A  (A6) 

which  is  in  the  form 

y  =  m  x  +  b  (A7) 

with 

m  =  1/A  (A8a) 

b  =  B/A  (A8b) 


From  Equations  (A5a)  and  (A5b)  it  can  be  seen  that  the 
solutions  for  a  given  by  Equation  (A7)  (if  any)  are  its  in¬ 
tercepts  with  the  unit  circle.  This  is  illustrated  in  Figure 
Al.  The  value  of  a  is  a  =  tan-1  (y/x),  where  {x,  y) 
are  the  coordinates  of  an  intercept. 

Two  special  cases  will  be  considered  first: 

•  m  =  0.  Then  y-b  and  x  =  ± (1  - b2)**  if  b  ^  1. 

•  b  =  0.  Then y-m  x  passes  through  the  origin,  so 

x2  +  y2  -  x2  +  (m  x)2  =  1 

x  =  ±1/(1  +  m2)'A,  y  =  m  x  -  ±m/(l  +  m2)^ 

In  both  of  these  cases  (if  b^l)  there  are  two  solutions. 

More  generally,  suppose-  m  ^  0,  b  &  0.  Draw  a 
perpendicular  from  the  origin  O  to  y  =  m  x  +  b,  in¬ 
tersecting  the  straight  line  at  P  =  (Xq,  y q).  The  perpen¬ 
dicular  line  has  the  equation 

y  =  (-  1/m)  x.  (A9) 

First,  let  us  find  the  distance  d  =  OP. 

m  x0  +  b  -  (-  1/m)  x0 

x0  =  -b  m/(m2  +  1)  (A  10a) 


y0  =  (-  1/m)  f-b  m/(m2  +  1)] 


=  b/(m2  +  1) 


(AlOb) 


y  ”  (-1/m)x 


y  -  mx  +•  s 


Figure  Al.  Equation  (A7)  displayed  on  the  unit  circle. 


d  =  (x02  +  y0 2)'A  =  [(b2  m2  +  b2)  /  (m2  +  1  )2]'A 


=  b/(m2  +  i)‘A 

There  are  three  possibilities: 


(A  10c) 


•  If  d>l,  then  Equation  (A7)  has  no  solution. 

•  If  d  =  1,  there  is  exactly  one  solution,  {Xq,  y0). 

•  If  d  <  1 ,  there  are  two  solutions. 


In  terms  of  A  and  B,  d  is  given  by 
d  =  B/(l  +  A2)'A 


(Alla) 


so  the  three  possibilities  above  can  be  expressed  in  terms  of 


/  (A,  B)  =  A2  +  1  -  B2. 


(Al  lb) 


Equation  (A7)  has  0,  1,  or  2  solutions,  depending  on 
whether /(A,  B)  is  negative,  zero,  or  positive,  respectively. 

When  f  (A,  B)>0  the  two  solutions  Rt  =  (xv  y,) 
and  R2  =  (x^  are  found  by  going  a  distance  s  = 
(1  -  d2)'A  from  point  P  in  both  directions  along  y  = 
m  x  +  b.  Equivalently,  they  can  be  found  by  applying 
the  identity  sin2 a  +  cos2 a  =  1  to  Equation  (A3).  The 
result  is 


a  -B  +  A(A2  +  1  -  B2)'A 
A2  +  1 


(A  12a), 


=  A  B  +  (A2  +  1  -  B2)‘A 
y>  A2  +  1 


(A  12b) 


lxij  l-B  +  A(A2+1-B2)»  J 


=  -B  -  A(A2  +  1  -  &)» 
A2  +  1 


=  A  B  ~  (A2  +  1  ~  B2)” 

yj  A2  +  1 


(A  12c) 


(A  13a) 


(A  13b) 


a  ^  =  tan 


A  B  -  (A2+l-B2)*  1 
-B  -  A(A2+1-B2)» J 
(A13c) 


The  occurrence  of/  (A,  B)  in  these  equations  should  be 
noted.  When  /(A,  B)  *  0  only  the  first  terms  in  both 
the  numerator  and  denominator  of  Equations  (A  12c)  and 
(A  13c)  remain,  a,  =  a^,  and  the  solution  for  possibili¬ 
ty  2  (discussed  after  Equation  (A  10c))  is  obtained. 

It  can  be  shown  that  at  and  a2  do  indeed  satisfy  Equa¬ 
tion  (Al).  Written  in  terms  of  A  and  B,  Equation  (Al) 
becomes 


A  =  cota.  +  B  csca  . 


(A  14) 

By  referring  to  a  unit  circle  diagram  it  can  be  seen  that 


=  0.  The  particular  situation  sinfi  =  0  occurs  when  the 
reflected  beam  is  in  the  zenith  direction.  Further  examina¬ 
tion  shows  that 


COtOtj 


5  .-B±A(A>+1-  B2)"  (A15>) 

>,  A  B  ±  (A2  +  l  -  B2)* 


If  v  =  0  and  <f>  +  ta  -  fi  =  0,  then  a  =  0°. 

If  v  =  0  and  <t>  +  ta  +  fi  =  0,  then  a  =  180°, 

If  v  =  180°,  then  a  =  180°  always. 


I  m  _ A2  +  1 _ 

>i  A  B  ±  (A2  +  1  -  B2)* 


(A  15b) 


The  cases  sinv  =  0  and  sinfi  =  0  will  be  discussed 
further  below. 


where  i  =  1  or  2  and  the  variable  signs  should  be  paired, 
+  with  +  and  -  with  -.  When  we  substitute  (A15a) 
and  (A  15b)  into  (A  14)  we  get 


A3. 3  sin/3  =  0 

B  is  defined.  Physically,  sinfi  =  0  means  there  is  no 
wave  slope,  so  clearly  a  has  no  meaning.  An  arbitrary 
value  such  as  0°  can  be  assigned  to  a  in  this  case. 


A  =  ~B±A(A2+ 1  -B2)‘A  +  B  (A2  +  1) 

A  B±(A2+1  -B2)'A  A  B±(A2+ 1  -B2)‘A 

=  A  [A  B  ±  (A2  +  1  -  B2)») 

A  B  ±  (A2  +  1  -  B2)'A 

=  A  .  (A16) 

A3.  Discussion 

There  are  several  special  cases  for  which  the  solution 
to  Equation  (Al)  cannot  be  performed  as  discussed  in  the 
preceding  section.  These  cases  will  be  taken  up  now. 

A3.1  cos</>  =  0 

B  =  0  so, 

A  sina  =  cosa, 

tana  =  1/A  =  1/cotv  =  tanv  .  (A17) 

Physically,  cos(j>  -  0  means  that  the  sun  is  at  the  zenith, 
so  the  direction  of  the  X2  and  y  axes  is  undefined.  Ar¬ 
bitrarily  we  may  pick  a  =  v  in  this  case,  but  it  is  not 
expected  to  occur  for  the  present  application. 

A3.2  sina  =  0 

Since  a  is  the  desired  solution  this  situation  is  not  ex¬ 
plicitly  apparent  in  advance.  When  sina  =  0  then  cota 
and  csca  are  undefined  (infinite)  and  Equation  (Al)  is  in¬ 
valid.  Returning  to  the  derivation  of  equation  (Al),  it  is 
seen  that 

sinv  sinfi  *  0.  (A  18) 

Physically,  it  is  clear  that  when  sina  ~  0  the  reflected 
beam  must  strike  the  image  plane  along  the  y-axis,  sinv 


A3. 4  cosco  =  0 

There  is  no  solution  for  a  from  the  original  geometric 
relations;  also,  B  is  undefined.  Physically,  costa  =  0  cor¬ 
responds  to  glancing  incidence.  An  arbitrary  value  such 
as  0°  can  be  assigned  to  a  in  this  case. 

A3.5  sin?  =  0 

A  is  undefined.  The  image  point  is  on  the  film  y-axis, 
which  can  only  happen  when  a  =  0°  or  a  =  180°. 
Refer  to  Section  A3.2. 

A3.6  sin#*  -0 

The  value  of  fi  does  not  appear  explicitly  in  Equation 
(Al),  so  presumably  the  solution  presented  in  Section  A2 
could  be  used.  But  it  was  shown  above  that  sinfi  -  0 
when  sina  =  0,  in  which  case  Equation  (Al)  is  not  valid. 
From  the  basic  geometrical  relationships,  when  sinfi  =  0, 

2  sina  sinfi  costa  -  0.  (A19) 

The  cases  sinfi  =  0  and  costa  =  0  have  already  been 
discussed.  These  possibilities  should  be  checked  for  first. 
If  sinfi  &  0  and  costa  &  0  when  sinfi  =  0,  then  sina 
=  0.  Physically,  sinfi  =  0  occurs  when  the  reflected  beam 
is  in  the  zenith  direction.  This  can  only  happen  when  a 
=  180°. 

Conclusion 

It  has  been  shown  that,  in  general,  there  are  two  solu¬ 
tions  for  a.  Both  solutions  are  mathematically  valid,  so 
the  choice  of  solution  in  a  particular  case  must  be  based 
on  further  analysis  of  the  physics  underlying  the  deriva¬ 
tion  of  Equation  (Al).  It  is  not  immediately  clear  how 
to  proceed. 
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Appendix  B:  Computer  program  improvements 


B1.  Introduction 

Some  changes  have  been  made  to  the  system  of  com¬ 
puter  programs  described  in  the  text  of  this  report,  for 
the  purpose  of  improving  the  results.  This  appendix 
describes  two  modifications.  The  first  of  them  is  a  pro¬ 
cedure  that  provides  a  better  approximation  to  the  original 
light  intensities  incident  on  the  film  than  simply  the 
digitized  pixel  values  themselves.  The  second  modifica¬ 
tion  is  a  revised  calculation  of  rms  slope  values  from  the 
Zx-Zy  histogram.  The  new  formulations  produce  some 
changes  both  in  the  programs  and  in  the  instructions  for 
use. 

B2.  Enhancements  to 
formulation 

B2.1  Correction  of  digitized  values  to 
obtain  original  luminance 

The  pixel  values  that  make  up  the  digitized  images  are 
related  to  the  reflected  brightness  field  that  illuminated 
the  film.  It  is  the  latter  (reflected  brightness)  that  should 
properly  be  used  to  construct  the  distributions  that  are 
analyzed  to  provide  information  on  the  sea  slope  distribu¬ 
tions.  The  former  (pixel  values)  is  what  is  actually  available. 
It  is  almost  certain  that  the  latter  values  are  not  directly 
proportional  to  the  former,  as  the  following  discussion 
shows. 

An  image  is  produced  when  light  of  intensity  IJx,  y) 
strikes  a  photographic  plate,  where  the  coordinates  x  and 
y  define  the  position  on  the  plate.  The  illumination  pro¬ 
duces  an  optical  density  distribution  D(x,  y)  when  the 
resulting  negative  is  developed  [10].  The  digitized  image 
is  produced  by  illuminating  the  negative  by  an  (ideally) 
uniform  light  intensity  /,.  The  intensity  of  the  light  that 
is  transmitted  is  given  by 

h  fa  y)  “  /;  lO-Wxy)  (Bl) 

In  the  NORDA  Remote  Sensing  Branch  IDSIPS  system 
the  intensity  disti  .i  ution  I2  (x,  yj  is  imaged  by  a  video 
camera,  whose  output  is  digitized  to  give  values  K  at  sam¬ 


pled  positions,  in  the  range  O^K  $$  255.  The  functional 
relationship  between  I2  and  K  is  unknown.  The  values 
AT,  which  describe  a  negative,  are  further  modified  to  give 
the  pixel  values  K  for  a  positive  image  by 

K  =  255  -  K.  (B2) 

To  perform  the  sun  glitter  analysis  it  is  necessary  to 
estimate  I0  from  K  .  Empirically,  it  is  found  that 

10- D  =  X  =  a  +  bK  +  cK?  (B3) 

describes  the  relationship  between  K  and  D  fairly  well. 
The  Hurter-Driffield  curve  that  describes  the  dependence 
of  D  on  exposure  for  photographic  film  ( D  -  log  E  curve) 
has  a  linear  region  given  by 

D  =  y  (log  I0t  -  log  i)  (B4) 

where  y  is  the  slope  of  the  linear  part,  t  is  the  exposure 
time,  and  i  is  a  constant  called  the  “inertia”  of  the  film 
[10]. 

If  we  assume  that  exposures  are  restricted  to  the  linear 
portion  of  the  D  -  log  E  curve,  we  can  estimate  I0  from 
K'  as  follows: 

•  Invert  Equation  (B2)  to  obtain  K. 

•  Use  Equation  (B3)  to  get  X  =  10~D  from  K. 

•  From  Equation  (B4), 

10-D  =  X  =  Al-  y  .  (B5) 

Therefore, 

l0  =  B/X!/y  .  (B6) 

The  choice  of  B  is  tantamount  to  the  choice  of  a  unit 
of  light  intensity,  so  B  can  be  chosen  for  convenience. 

The  only  unknowns  in  the  above  procedure  are  the  three 
coefficients  in  Equation  (B3),  a,  b.  and  c.  Table  Bl  lists 
a  set  of  calibration  values  that  were  used  to  obtain  the 
coefficients.  A  least-squares  fit  to  the  data  of  Table  Bl  gave 

a  =  0.10138  x  10- 1 

b  =  0.97293  x  10~) 

c  =  0.21483  x  10-4 
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Table  B1 .  Step  wedge  calibration  values. 


Step 

0 

10-D 

K 

1 

0 

1 

186 

2 

0.1 

0.794 

172 

3 

0.2 

0.631 

159 

4 

0.38 

0.417 

117 

5 

0.59 

0.257 

81 

6 

0  83 

0.148 

53 

7 

1.04 

0.0912 

39 

8 

1.28 

0.0525 

29 

9 

1.50 

0  0316 

21 

10 

2.27 

0.00537 

12 

D  = 

Density  ot  step 

K  = 

Resulting  digitized  value 

B2.2  RMS  slope  calculation 

In  the  computer  program  described  in  the  text  of  this 
report,  rms  slope  values  were  calculated  in  the  following 
way:  Two  “slices”  of  the  distribution  of  brightness  over 
components  of  slope  were  taken,  one  along  the 
upwind  =  downwind  axis  and  the  other  along  the  cross- 
wind  axis.  Each  slice  was  taken  to  be  a  separate  univariate 
distribution,  from  which  the  rms  slopes  (upwind  and  cross- 
wind)  were  calculated. 

To  state  the  above  formally,  let  the  bivariate  distri¬ 
bution  be  called  p  (x,  y).  Since  p  (x,  y)  is  a  probability 
density, 

j  J  P  (x,  y)  dx  dy  =  1  (B7) 

where  the  integral  is  over  the  appropriate  region  of  R2- 
(There  is  an  equivalent  formulation  for  discrete  distribu¬ 
tions.)  The  two  slices  are  Ap(x,  0)  and  B  p(0 ,  y)  where 
A  and  B  are  chosen  to  normalize  the  integral  over  x  or 
y ,  respectively.  The  rms  values  were  calculated  as  the 
square  roots  of 

(B8) 


and  the  corresponding  integral  over  the  other  slice.  In 
Equation  (B8)  m  is  the  mean  value  of  x,  (In  fact,  m  = 
0  was  assumed  in  SUBROUTINE  RMSARY.) 

From  usual  statistical  theory,  it  might  be  assumed  that 
Equation  (B8)  should  be  replaced  by 

\  \  (x  -  m)2  p  ( x ,  y)  dx  dy  (B9) 

(and  the  equivalent  for  the  y  component).  The  m  in  Equa¬ 
tion  (B9)  is  not  necessarily  the  same  as  the  m  in  Equation 
(B8).  The  discussion  by  Cox  and  Munk  does  not  make 
it  clear  which  formulation  is  correct  [1,  2],  Furthermore, 
for  some  distributions  (such  as  the  bivariate  normal 
distribution  with  no  correlation  between  x  and  y)  the 
results  are  no  different.  However,  for  completeness  the 
formulation  of  Equation  (B9)  was  implemented  and  tried. 

B3.  Changes  to  programs  and 
user  instructions 

B3.1  Original  intensity  calculations 

Equations  (B2),  (B3),  and  (B6)  were  implemented  in  a 
new  subroutine,  INTENS.  The  main  program  was 
modified  to  a  new  version,  GMAIN3,  which  calls  INTENS 
and  contains  the  other  necessary  changes.  INTENS  was 
compiled  into  the  same  RL,  GSUBSRL,  as  the  other 
subroutines.  GMAIN3  was  compiled  into  GMAIN3RB. 
The  prepared  program  has  the  file  name  GLITTER3. 

GLITTER3  requires  the  same  :FILE  statements  as 
GLITTER2.  The  only  change  is  that  one  extra  input  item, 
7,  is  needed.  It  goes  in  a  separate  record  immediately 
following  focal  length  and  picture  height. 

B3.2  Revised  RMS  slope  calculation 

The  changes  for  the  calculation  described  in  Section 
B2.2  are  in  a  new  version  of  the  main  program,  GMAIN4. 
GMAIN4  is  a  revision  of  GMAIN3.  Like  the  latter  it 
calls  INTENS,  and  7  is  a  required  input.  GMAIN4  was 
compiled  into  GMAIN4RB.  The  prepared  program  has 
the  file  name  GLITTER4.  (RMSARY  is  no  longer  used.) 
The  instructions  for  use  are  exactly  the  same  as  for 
GLITTER3. 
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